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Abstract 

A numerical model is developed to examine laminar 
flame spread and extinction over a thin solid fuel in low- 
speed concurrent flows. The model provides a more precise 
fluid— mechanical description of the flame by incorporating 
an elliptic treatment of the upstream flame stabilization 
zone near the fuel burnout point. Parabolic equations are 
used to treat the downstream flame, which has a higher flow 
Reynolds number . The parabolic and elliptic regions are 
coupled smoothly by an appropriate matching of boundary 
conditions. The solid phase consists of an energy equation 
with surface radiative loss and a surface pyrolysis 
relation. Steady spread with constant flame and pyrolysis 
lengths is found possible for thin fuels and this facili- 
tates the adoption of a moving coordinate system attached 
to the flame with the flame spread rate being an eigen- 
value. Calculations are performed in purely forced flow in 
a range of velocities which are lower than those induced in 
a normal gravity buoyant environment. Both quenching and 
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blowoff extinction are observed. The results show that as 
flow velocity or oxygen percentage is reduced, the flame 
spread rate, the pyrolysis length, and the flame length all 
decrease, as expected. The flame standoff distance from 
the solid and the reaction zone thickness, however, first 
increase with decreasing flow velocity, but eventually 
decrease very near the quenching extinction limit. The 
short, diffuse flames observed at low flow velocities and 
oxygen levels are consistent with available experimental 
data. The maximum flame temperature decreases slowly at 
first as flow velocity is reduced, then falls more steeply 
close to the quenching extinction limit. Low velocity 
quenching occurs as a result of heat loss. At low veloci- 
ties, surface radiative loss becomes a significant fraction 
of the total combustion heat release. In addition, the 
shorter flame length causes an increase in the fraction of 
conduction downstream compared to conduction to the fuel. 
These heat losses lead to lower flame temperatures, and 
ultimately, extinction. This extinction mechanism differs 
from that of blowoff, where the flame is unable to be 
stabilized due to the high flow velocity. 


IX 



Table of Contents 


Abstract 


Table of Contents 


in 


List of Figures 

Nomenclature 

Chapter 1. Introduction 

1. Classes of Solid Combustion 

2. Solid Phase Combustion Modelling Techniques 

3. Overview of Early Relevant Solid 
Combustion Models 

Microgravity Combustion of Solids: 

Modelling and Experiments 

1.5. Combustion Experiments in Space 

1.6. Summary 

Chapter 2. Theoretical Formulation 


1 

1 

1 


1.4 


2.1. Elliptic Region 

2.1.1. Flow Field Equations 

2.1.2. Flow Field Boundary Conditions 

2.1.3. Nondimensional Equations and 

Boundary Conditions 

2.1.4. Energy Equation 

2.1.5. Energy Equation Boundary Conditions 

2.1.6. Nondimensional Energy Equation 

and Boundary Conditions 

2.1.7. Fuel and Oxygen Species Conservation 

Equations 

2.1.8. Fuel and Oxygen Species Boundary Conditions 

2.1.9. Nondimensional Species Equations 

and Boundary Conditions 

2.2. Parabolic Region 

2.2.1. Nondimensional Parabolic Region Equations 

2.2.2. Nondimensional Parabolic Region 

Boundary Conditions 

2.3. Coupling Between Elliptic and Parabolic Regions 

2.4. Solid Phase Equations 


vm 


2 

7 

8 

13 

18 

20 

22 

25 

25 

27 

28 

31 

32 

33 

34 

34 

35 

36 
36 

36 

37 

38 


iii 



Chapter 3 . Numerical Model 44 

3.1. Discretization of Convective Terms 45 

3.2. Treatment of Source Terms 47 

3.3. Staggered Grid System 48 

3.4. The SIMPLER Algorithm 50 

3.5. Domain Structure 52 

3.6. Parabolic Region 54 

3.7. Computer Usage 55 

3.8. Property Values 56 

Chapter 4 . Results 57 

4.1. Detailed Flame Characteristics 

For One Case 58 

4.1.1. Gas Phase Profiles; 5 cm/s, 15% 0 2 58 

4.1.2. Solid Phase; 5 cm/s, 15% 0 2 70 

4.2. Comparison With Experiment 73 

4.3. Parametric Comparison of Theoretical Results 78 

4.3.1. Global Results: Spread Rates and 

Maximum Temperature 78 

4.3.2. Parametric Comparison of Flame Structure 81 

4.3.3. Flame Stand-off Distance and Thickness 83 

4.3.4. Solid Phase Parametric Results 89 

4.3.5. Definition of Various Length Scales 91 

4.4. Extinction Boundary 93 

4.4.1. Differences Between Blowoff and Quenching 97 

4.4.2. Incomplete Combustion 104 

Chapter 5. Conclusions 108 

References 

Appendix 


IV 



Fig. 1 
Fig. 2 


List of Figures 


Flame spread over solids, (p.4) 

Schematic of problem showing elliptic and para- 
bolic regions, (p.23) 

Fig. 3 Elliptic/parabolic matching boundary conditions. 
(p.39) 

Fig. 4 Geometry of solid fuel, (p.40) 

Fig. 5 Grid structure for main variables Q, and u and v- 
velocities. (p.49) 

Fig. 6 Variable grid structure used in this problem, 

(p • 53) 

Fig. 7 Computational test matrix, (p.59) 

Fig. 8 Nondimensional isotherms. Uj„=5 cm/s, X 0 =15% 

(p.60) 

Fig. 9 Concentration profiles. U^=5 cm/s, X 0 =15% (p.62) 

Fig. 10 Various flame length depictions. 1^=5 cm/s, 

X 0 =15% (p . 64 ) 

Fig. 11 Velocity vector plot. "0^=5 cm/s, X 0 =15% (p.66) 

Fig. 12 Velocity vector plot (lab. coordinates) . "0^=5 

cm/s, X 0 =15% (p.68) 

Fig. 13 Streamlines in_ flame (top) and lab. (bot . ) 

coordinates. U^,=5 cm/s, X 0 =15% (p.69) 

Fig. 14 Isobars. "0^=5 cm/s, X 0 =15% (p.71) 

Fig. 15 Solid phase profiles: incident heat flux, surface 
temperature, fuel thickness, and blowing veloci- 
ty. 0^=5 cm/s, X 0 =15% (p.72) 

Fig. 16 Comparison of calculated fuel reactivity levels 

(a) with experimental flame shape (b) . U„~5 

cm/s, X 0 =15% (p.75) 

Fig. 17 Dimensional profiles. "0^=5 cm/s, X 0 =15% (p.76) 


v 



Fig. 18 Comparison of model and experiment at a free 
stream velocity of approximately 5 cm/s. (p.77) 

Fig. 19 Spread rate and T max as a function of relative 
velocity, (p.79) 

Fig. 20 Fuel reactivity contours at 15% 0 2 . Weakest 

(outermost) contour is 10~ 6 g/cm 3 /s. A factor of 
10 separates contours, (p.82) 

Fig. 21 Dimensional fuel reactivity contours at 15% 0 2 . 

Weakest (outermost) contour is 10 -6 g/cm 3 /s. A 
factor of 10 separates contours, (p.84) 

Fig. 22 Fuel reactivity contours at 5 cm/s. Weakest 

(outermost) contour is 10 -6 g/cm 3 /s. A factor of 
10 separates contours, (p.85) 

Fig. 23 Dimensional fuel reactivity contours at 5 cm/s. 

Weakest (outermost) contour is 10 -6 g/cm 3 /s. A 
factor of 10 separates contours, (p.86) 

Fig. 24 Stand-off distance, l s , and flame thickness, l t . 
(p. 88) 

Fig. 25 Solid temperature and heat flux profiles at 15% 
0 2 . (p . 90) 

Fig. 26 Solid temperature and heat flux profiles at 5 
cm/s. (p.92) 

Fig. 27 Flame (1 F ), pyrolysis (l p ) , and preheat (l h ) 
lengths, (p.94) 

Fig. 28 Extinction boundary, (p.96) 

Fig. 29 Extinction boundary plotted in regular coordi- 
nates. (p.98) 

Fig. 30 Blowoff extinction sequence. Fuel reactivity 
contours are shown. The outermost (smallest) 
contour is 10 -7 g/cm 3 /s. A factor of ten sepa- 
rates contours, (p.99) 

Fig. 31 Isotherms for a near quench point (a) and near 
blowoff point (b) at 15% 0 2 . (p.101) 

Fig. 32 Variation of flame spread rate and maximum 
temperature along the extinction boundary. 
(p.103) 


vi 


Fig. 33 Flame structure effects near blowoff. Dotted 
lines are fuel mass fraction. Solid lines are 
fuel reactivity contours, (p.105) 

Fig. 34 Fuel flux in x-direction at different free stream 
velocities, (p.106) 


Vll 



Hi h| h w HI ,£) ooo|ti h| h 3- t-< f fl tr o >n| h> n h| n n| D| o o o| o| w| >| 

*8 •< h* 8 co h wWiQ^H**P- 


Nomenclature 


Solid Phase Pre -Exponential Factor (3.8xl0 7 cm/s ) 

Gas Phase Pre-Exponential Factor (l.OxlO 13 cm 3 /g/s) 
Gas Phase Specific Heat (0.33 cal/g/K) 

Solid Phase Specific Heat (0.30 cal/g/K) 

Damkohler Number, a*f>*B g /U R 

Diffusion Coefficient of i (i=F or 0), 0^0* = (T/T*) 2 
Reference Diffusion Coeffient (2.13 cm 2 /s) 

Gas Phase Activation Energy (2.7xl0 4 cal/gmol) 

Eg / R / (45.3) 

Solid Phase Activation Energy (3.0xl0 4 cal/gmol) 

Es / R / T m (50.3) 

Oxygen/ Fuel Stoichiometric Mass Ratio (1.185) 

Gravity Level (0.0 cm/s 2 ) 

Grashof Number, a*g/U R (0.0) 

Solid Thickness, h/x R 

Latent Heat of Fuel (-180 cal/g) 

Z/C s /T eo (- 2 . 00 ) 

L Lewis Number of i ( i = F or 0) , a* /D * (1.0) 

Mass Flux from Fuel, iii/p”VU R 
Pressure, (P-F.J /p*/U R 
Ambient Pressure (1 atm) 

Prandtl Number, v*/a* (0.7) 

Heat of Combustion (4.00xl0 3 cal/g) 

Q/Cp/lL (40.4) 

Energy Source Term, -Q© F 
Heat Flux to Solid, K(0T/8y) w 

Ideal Gas Constant (1.989 cal/gmol/K) 

Distance Along Fuel Surface, "s/xr 

Temperature, T/T^ 

Ambient Temperature (300 K) 

Reference Temperature (1250 K) 

T* T*/ T,*, (4.167 ) 


viii 


T f Adiabatic Flame Temperature (2200 K) 

T l Temperature at which L is given ( 300 K ) 

T l Tl/T^ (1.0) 

T s Solid Temperature 
T s 

u Velocity in x-direction, u/U R 
U B Reference Buoyant Velocity, (a*g) 1/3 
U R Reference Velocity, u^ + Ub-Vj, 

Forced Flow Velocity 
v Velocity in y-direction, v/U R 
V F Flame Spread Rate 
x x-coordinate, x/x R 
x R Reference Length, a*/U R 
y y-coordinate, y/x R 
Y f Fuel Mass Fraction 
Y 0 Oxygen Mass Fraction 
Y 0/ o» Ambient Oxygen Mass Fraction 

a* Reference Thermal Diffusivity (2.13 cm 2 /s) 
r Solid Phase Radiation Parameter, a£T^/ (p’*C^ p u R ) 

£ Solid Emissivity (1.0) 

K Gas Thermal Conductivity, k/k* = tVt"* 

k* Reference Thermal Conductivity (1.93X10" 4 cal/cm/s/K) 
P Gas Viscosity, (I/jI*=T/T* 

)i* Reference Gas Viscosity ( 4 . 10 x 10“ 4 g/ cm/s ) 

H Solid Phase Parameter, p s V F C s / p*/ U R /C P 

p Gas Density, p/p - * 

p* Reference Density ( 2 . 75 x 10" 4 g / cm 3 ) 

P<» Ambient Density ( 1 . 15xl0 -3 g/cm 3 ) 

Poo Density, p^/p* (4.167) 

Ps Solid Density (0.263 g/cm 3 ) 

o Stefan-Boltzmann Constant ( 1 . 356 xlO" 12 cal / cm 2 / K 4 / s ) 
x Fuel Half-Thickness (0.0038cm) 

% Fuel Source Term, -Dap 2 Y F Y 0 exp ( -E g /T ) 

Oxygen Source Term, fcbp- 


IX 



Subscripts 


B Buoyant 

f Flame 

F Fuel or Flame 

g Gas Phase 

h Heat up 

i Species i = F or 0 

L Refers to Latent Heat 

max Maximum 

min Minimum 

0 Oxygen 

p Pyrolysis, Pressure 

par Refers to parabolic region 
prod Products 
R Reference 

s Solid Phase 

S Stand-off 

T Thickness 

w Value at Fuel Surface 

y Along the y-direction 

oo Value at Far Field 


Superscripts 

* Evaluated at T* 

Note: A bar above a variable indicates that it 

dimensional quantity. 


is a 


x 



Chapter 1 . Introduction 

The study of combustion in a microgravity environment 
offers a chance to examine aspects of flames without the 
complicating influence of natural convection . 1 Recently, 
interest in microgravity combustion has been steadily in- 
creasing. This is due in large part to experiments already 
performed in microgravity, each of which has revealed new 
and exciting combustion phenomena . 2 

While the body of scientific knowledge of microgravity 
combustion has been getting larger, there is still an ever 
increasing series of new questions which arise. Therefore, 
the demand for microgravity facilities (such as the Space 
Shuttle, drop towers, and Keplerian trajectory aircraft) 
has increased considerably. 

In addition to the scientific reasons, a practical 
concern of studying combustion in microgravity is space- 
craft fire safety. In a recent paper detailing past and 
present fire safety practices aboard spacecraft, it is 
concluded, among other things, that "the microgravity 
research community has much to offer advanced spacecraft 
fire safety ." 3 

As is true with any scientific research effort, an 
experiment should be complemented by a theoretical investi- 
gation. The theoretical side of the research is often 
guided by the experiment, but likewise can aid in the 
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selection of the experiment. When seen in this light, the 
theoretical effort is simply the formalization and expan- 
sion of the knowledge obtained in the experiment. 

Additionally, a sound theory is useful in predicting 
results when an experiment is infeasible. This is espe- 
cially evident in microgravity studies. Any earthbound 
facility is limited in its ability to provide a micro- 
gravity environment, and of course a space experiment is 
very costly. Thus, a capable theory becomes a valuable 
tool both in understanding and predicting the physics of 
the problem. 

1.1. Classes of Solid Combustion 

Solid combustion is a broad field. Its study is made 
tractable by dividing it into several classes. First, 
almost all common fires are composed of diffusion flames. 
A diffusion flame initially separates the fuel and oxidiz- 
er, which diffuse into the flame zone and react, hence the 
name . * 

Flames can be stationary or spreading. A stationary 
flame is established when only a fixed area is exposed to 
the flame. On the other hand, a spreading flame, which is 
of more practical concern, greatly depends on the rate at 


*In a premixed flame, on the other hand, the fuel and 
oxidizer are initially mixed. The combustion of a double- 
base solid propellant is an example. 
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which the solid fuel can be heated up to its pyrolysis 
temperature. Thus, any prediction of a spreading flame 
needs to consider the solid phase. 

A further division is made by considering the direc- 
tion of the oxidizer flow relative to the direction of 
flame spread, shown in fig. 1. In opposed flow flame 
spread, the flame spreads in the direction opposite the 
flow (fig. la), while in concurrent flow flame spread, the 
flame spreads in the same direction as the flow (fig. lb) . 
The oxidizer stream can consist of forced flow, buoyant 
flow, or a combination of the two. On earth, a hot flame 
generates buoyant flow up. Thus, a flame spreading 
downward spreads in opposed flow mode, and a flame spread- 
ing upward spreads in concurrent flow mode. 

Concurrent flow flame spread has received much atten- 
tion, since fires of practical interest are most hazardous 
in this configuration. Many experiments and analytic 
studies have been performed on concurrent flow flame 
spread. Detailed numerical predictions have been somewhat 
hindered by the complicated nature of the problem, since 
the flames are often large, turbulent, and radiant. 

Whether opposed flow or concurrent flow, certain 
fundamental mechanisms apply. The flame transfers heat to 
the solid by conduction and radiation. This heat causes 
the solid to gasify, and the gaseous fuel then flows toward 
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the flame. In the flame zone, fuel and oxygen come 
together, react, and sustain the flame. The flame is 
stabilized in a position that balances the complicated 
flows of heat and mass. 

To get a complete understanding of these flames, we 
are interested not only in their structure and spread 
characteristics, but also in extinction (or how they go 
out) . Everyone is familiar to some degree with the effect 
of flow velocity on flames. For example, gently blowing on 
a weak fire usually intensifies it. However, blow too 
hard, and the fire will go out . This "blowoff " is due to 
the fact that too much oxygen is supplied to the flame. 
The reaction rate in the flame zone simply cannot keep pace 
with the high flow rate of oxygen, and the flame is 
literally blown away. Blowoff extinction has been studied 
extensively . 

On the other side of the coin, what would happen to a 
flame if the flow of oxygen were reduced steadily, perhaps 
all the way to zero? In the gravitational field of earth, 
we cannot answer this question, since buoyant flows (on the 
order of tens of centimeters per second) will always 
accompany any flame. By studying flames in microgravity, 
we can endeavor to answer this and other questions. Will 
the flame burn continuously, or will it eventually go out? 
Why? What is the effect of low-speed flow on the flame? 
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A theoretical model 4 suggests that the flame will eventual- 
ly go out as the flow velocity is reduced. The reason is 
that heat loss becomes important. The heat loss reduces 
the flame temperature below the point that combustion can 
be sustained. 

In this effort, concurrent-flow flame spread over a 
thin fuel in zero gravity* is to be studied. The flow 
field is generated by purely forced convection. (There is 
no buoyant flow since gravity is assumed to be zero.) 
Whenever we say "thin fuel" we mean hydrodynamically and 
thermally thin. Hydrodynamically thin means that the 
thickness of the fuel is always much smaller than the 
distance of the flame from the fuel surface (or, equiva- 
lently, the gas phase length scale) . Thermally thin means 
that the rate of conduction in depth of the solid fuel is 
much faster than the heat up rate of the gas phase along 
the fuel surface (effectively, this means that the tempera- 
ture is constant across the fuel thickness) . 

Before proceeding with a detailed explanation of the 
problem, a brief presentation of earlier related work on 
solid phase flame spread is given below. The following 
discussion is not intended to be an exhaustive review, but 

In a spacecraft, a microgravity environment exists. 
We consider strictly zero gravity, however, as a starting 
point. To model microgravity effects, buoyant and inertial 
forces should be considered. 


7 


instead outlines the historical framework which provided 
the motivation for this work. 

1.2. Solid Phase Combustion Modelling Techniques 

The earliest models relied on heat transfer arguments 
to obtain expressions for global parameters (e.g., spread 
rate) in terms of other known quantities. Despite their 
simplicity, trends were predicted fairly well. (Today, 
more sophisticated analytical predictions of piloted 
ignition and flame spread using fuel and environment 
parameters as inputs are used to try to cover a wide range 
of conditions. 5 ) 

More sophisticated models solved conservation equa- 
tions to obtain a better understanding of the flames. The 
earliest of these solved parabolic (or boundary layer type) 
equations. Sometimes, similarity solutions existed. 
Otherwise, these equations could be solved by a numerical 
marching scheme, i.e., given an initial profile, the 
downstream profile at the next grid point could be found 
immediately. Furthermore, the "flame sheet" approximation 
was made. This assumed that as soon as the fuel and 
oxidizer met, they reacted (infinitely fast) . This 
neglected the detailed chemical reaction effects. 

Later models began to include real, finite-rate 
chemical reaction effects . The added difficulty was in the 
evaluation of highly non-linear Arrhenius rate relations. 
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Finally, models included the complete (elliptic) 
conservation equations, together with finite-rate reaction 
expressions. The solution of elliptic equations is 
considerably more difficult than parabolic equations, since 
any point in the domain is, in general, affected by all 
surrounding points, not only upstream points. Early models 
of this type assumed a velocity field, but eventually, even 
the full Navier-Stokes equations were solved. 

Some studies of flames above solid fuels were simpli- 
fied by considering the case where the flame was station- 
ary. How is a stationary flame established? The trick is 
to allow only a small area of fuel to be exposed, and 
assume the rate of fuel regression (or burnout) is small. 
A good simulation of such a flame is made by considering a 
flame over a fuel soaked wick, or a flame above a porous 
burner through which fuel is forced. In either case, of 
course, no burnout can occur. Stationary flame studies are 
well-suited to examining the flame structure as a function 
of environmental parameters, but cannot examine flame 
growth and spread. In the following discussion, both 
stationary flame and flame spread models will be presented. 
1.3. Overview of Early Relevant Solid Combustion Models 

An approximately chronological summary of combustion 
modelling research relevant to this work is presented. As 
was mentioned previously, only a cursory overview is given. 
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The reader is advised to examine several review arti- 
cles 6 ' 7 for a more complete early history, as well as a 
more recent review . 6 

An early work 9 examined the stationary flame structure 
of a fuel soaked wick both with an experiment and a 
numerical model. In the model, laminar, boundary layer 
type equations were solved to obtain the flame structure in 
the flame and thermal plume. The initial profile for the 
equations was provided by a similarity solution which 
exists over the region where pyrolysis was occurring (i.e., 
above the wick) . Infinitely fast kinetics (flame sheet 
approximation) were assumed. While the prediction approxi- 
mated the experimental results fairly accurately, the model 
was simple by today's standards. Yet, by the authors' own 
admission, it was too awkward to use for direct flame 
spread analysis, given their computational capabilities. 

A predictive model of flame spread in concurrent flow 
in several different geometries (floor, ceiling, and wall, 
e.g.) was produced . 10 The solid phase was unsteady, and the 
gas phase quasi-steady, i.e., the gas phase response time 
was much shorter than the characteristic solid time. The 
gas phase equations were two-dimensional, laminar, and 
parabolic (boundary layer type) . Radiation was neglected, 
and a flame sheet approximation was made. Using various 
assumptions, the solid and gas phase equations were 
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decoupled and a similarity solution was found for all 
equations except for the species equation downstream of the 
pyrolysis zone (a similarity solution for the ceiling 
configuration was not found, however) . 

One of the limitations of these works was the 
neglection of finite-rate chemical kinetic effects. Since 
a flame is a complicated reacting system, the inclusion of 
an improved chemical model is necessary. The effects of 
temperature, fuel and oxidizer concentration, flow veloci- 
ty, and mass diffusion then can be seen to influence flame 
spread and extinction characteristics through the flame 
chemistry. Including finite-rate chemistry, the elliptic 
two-dimensional, quasi-steady gas phase energy and species 
equations were solved . 11 The flow field was assumed uni- 
form, thus eliminating the need to solve the complicated 
Navier-Stokes equations. The solid fuel equation neglected 
conduction ahead of the flame, but included an Arrhenius 
type pyrolysis relation and unsteady st>lid phase heat up. 
Opposed flow flame spread over a thin solid fuel was 
examined. By considering opposed flow flame spread, the 
problem is immediately simplified because the flame is not 
only shorter*, but the spread process only depends on the 
flame structure at the anchor point, whereas in concurrent 

*The shorter flame can make the effects of radiation 
and turbulence less important. 
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flow, the entire flame and thermal plume should be consid- 
ered to predict spread rate. Flame spread and blowoff 
extinction were obtained, in good agreement with experi- 
ment. This model provided the basis for many of the more 
recent ones, including this work, which are able to utilize 
improved computer capability and performance. 


Another attempt to go beyond the flame sheet assump- 
tion was made . 12 A stationary flame calculation, which 
solved the laminar boundary layer equations for velocity, 
temperature, and species, included a finite— rate chemical 
reaction term. As done in earlier work, the initial 
profile was provided by the similarity solution which 
exists in the pyrolysis zone. The finite rate chemistry 
effects led to a shorter flame length and fuel preheat 
distance. In addition, there was the possibility of fuel 
vapor escaping from the flame zone unreacted (escaped 
pyrolyzates or incomplete combustion) , which could be 
quantified . 


Modelling continued to move toward a more thorough 
formulation. The complete elliptic set of equations, 
including the Navier-Stokes equations, were used to examine 
the mixed flow* combustion of a vertical fuel plate imbed- 
ded in an inert substrate . 13 Laminar, steady flow was 

Mixed flow includes any combination of buoyant and 
forced flow. 
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assumed. The fuel plate was 4 centimeters long, and the 
inert substrate extended 1 centimeter in front of the fuel. 
Fuel burnout was assumed to be negligible, so a stationary 
flame existed. The small fixed length of the fuel permit- 
ted computation, since a fuel of longer extent would 
require a prohibitively large computational time. A simple 
energy balance was used to model the rate of production of 
pyrolysis products from the fuel. The equations were 
solved in dimensional form, and results demonstrated the 
importance of the pressure field in controlling the flow 

near the flame stabilization point. 

A similar examination of a stationary flame stabilized 
at the leading edge of a fuel plate was carried out. 
The main difference from the previous work was the fact 
that the fuel began immediately (i.e., there was no inert 
plate extension) and the fuel plate extended all the way 
downstream. Again, the full two-dimensional laminar 
elliptic equations were solved (here, in nondimensional 
form) including finite-rate chemical reaction. The effect 
of flow velocity on the structure of the flame was examined 
(buoyant flow was neglected) . At first, the flame envel- 
oped the fuel completely, but as the flow velocity was 
increased, the flame retreated downstream, and was eventu- 
ally blown out of the computational domain. The importance 
of using elliptic equations to model the stabilization of 
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the flame was demonstrated. Furthermore, at arbitrarily 
low flow velocities, a flame which was dimensionally small 
but kinetically strong was found to exist. This was 
because, with the neglection of any heat loss mechanism, 
the flame had plenty of time for the combustion reactions 
to proceed, but merely reduced its size. 

A concurrent flow flame spread model examined the 
combustion of a thin fuel. 1 ^ This laminar, unsteady 
formulation was able to predict the initial transitory 
flame growth period followed by steady flame propagation. 
The laminar, two-dimensional, elliptic energy and species 
conservation equations were solved numerically. A simple 
finite-rate chemical reaction was assumed. Solution of the 
Navier— Stokes equations was avoided by prescribing a 
velocity field, in this case, a constant property Hagen- 
Poiseuille flow. Relatively high velocity forced flows 
(i.e., 60 cm/s and up) were studied. Some agreement with 
experimental data was obtained, all the way up to the 
blowoff limit. 

1.4. Microgravity Combustion of Solids: 

Modelling and Experiments 

Up to this point, the discussion has presented work 
which was not necessarily interested in the effects of 
microgravity. If gravity was neglected in a work, a high 
forced flow velocity was substituted. In ref. 14, there 
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were hints of the utility of studying that type of combus- 
tion in a reduced gravity environment. At low forced flow 
velocities (sub-buoyant) , the flame was very small, and 
could be more susceptible to heat loss effects. Thus, the 
solid fuel was allowed to lose heat through black body 
radiation . 16 Now, as the flow velocity was decreased, the 
flame eventually went out due to the increased relative 
importance of heat loss. Thus, both quenching and blowoff 
extinction were observed. 

At this time, the importance of an elliptic treatment 
of the flame stabilization zone was recognized. This is 
due to the fact that in the stabilization zone, the thermal 
length is the appropriate scale (meaning the product of the 
Reynolds and Prandtl numbers is unity, or Re x Pr=l) , and 
thus the Reynolds number in the flame stabilization zone 
was order unity* . When the Reynolds number is order unity , 
diffusion of mass, momentum, and heat in both the stream- 
wise and cross-stream directions is important. Thus, 
elliptic type equations result. A unified presentation of 
quenching and blowoff extinction for several flame systems 
was made . 


*The Prandtl number is order unity for most common 
gases. It is a measure of the rate of diffusion of 
momentum to the rate of conduction of heat for a given 
fluid . 
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In more concrete terms, the thermal length is found by 
considering a simple convection/conduction balance. As the 
cold oxygen stream flows into the flame and is warmed to 
the flame temperature, its rate of change of energy is 
given by the expression puC p (Tf-T^) , where p, u, and C p are 
the gas density, velocity and specific heat respectively. 
T f is the flame temperature, and T^ is the ambient tempera- 
ture. This flow is heated by conduction from the flame, 
given approximately as K(T f -T oc )/x, where k is the thermal 
conductivity of the gas and x is thermal length. Equating 
the two expressions yields the thermal length, x=k/ (puC p ) . 
Then, defining a Reynolds number based on x, this last 
expression is equivalent to Re x Pr=l . 

The experimental efforts examining flame spread over 
solid fuels in microgravity began to expand. Quiescent 
flame spread over a thin fuel in microgravity* was stud- 
ied . 18 Experiments were conducted at atmospheric pressure 
over a range of oxygen percentages from pure oxygen down to 
the limiting value. Among the findings, quenching extinc- 
tion was observed in microgravity. It was found to be 
quite different from blowoff extinction encountered in 
normal gravity. The flame structure and low flow veloci- 

*A flame spreading in a quiescent, microgravity 
environment spreads in opposed flow mode. This is clear in 
a flame-fixed coordinate system, where the oxidizer feeds 
into the flame in the opposed flow configuration. 
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ties verify the need for an elliptic system of equations to 
model the problem in that the Reynolds number and Peclet 
number are order unity. This is especially true for slow 
flames in microgravity, where the flame stabilization 
region is a large percentage of the overall flame zone. 

The previous results were combined with a study 
examining the effect of low speed flow on these flames 
spreading in microgravity 19 to clarify the role of convec- 
tion on opposed flow flame spread and extinction . 20 An 
extinction boundary was generated. Another more recent 
investigation 21 filled in additional data. For a flame 
burning at a given oxygen percentage, as the characteristic 
velocity was increased, the flame would eventually be 
blown off and as the characteristic velocity was decreased, 
the flame eventually was quenched. Between the two limits, 
a point where the flame would spread fastest . 
This point could be at a velocity below that due to normal 
gravity buoyant flow. Hence, the most hazardous condition 
from a fire safety standpoint could happen if buoyancy is 
completely removed and a small forced flow is applied. In 

The characteristic velocity represents the rate at 
which oxygen is convected into the flame zone. It is best 
visualized in flame-fixed coordinates. For example, for 
opposed flow flame spread, the characteristic 
velocity is the flame spread rate. For a downward burning 
flame in a gravitational environment, the characteristic 
velocity is the magnitude of the buoyant flow at the flame 
stabilization point, plus the spread rate. 
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a spacecraft, therefore, ventilation currents may establish 
this scenario. The importance of further work was clear. 

On the modelling front, a sophisticated prediction of 
opposed flow flame spread characteristics at microgravity 
was made . 22 Imposed flow velocities from zero (no flow, 
quiescent spread) all the way up to the blowoff limit, were 
studied. This steady, laminar, two-dimensional model nu- 
merically solved the full Navier-Stokes equations together 
with elliptic energy and species concentration equations. 
Radiation from the gas and solid phase, and their interac- 
tion, was modelled. Sample flame structures and radiation 
profiles were given. It was found that flames in micro- 
gravity are radiatively controlled, a phenomenon that would 
be masked in normal gravity. Specifically, results showed 
that including gas phase radiation greatly changed the 
flame structure in that a cooler, smaller flame results. 
Although, by including only surface radiation (and neglect 
ing gas phase radiation completely) flame spread rates and 
extinction trends were largely unchanged. This suggests 
that the important aspects of low-speed flame spread 
modelling can be captured without including the very 
difficult gas phase radiation treatment, but instead only 
the essential heat loss mechanism given simply by solid 
phase radiation. However, it is clear that gas phase 
radiation is necessary to complete the task of quantita- 
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tively predicting flames, when detailed chemistry and soot 
production need to be included. 

Very recently, the last model has been upgraded to 
examine unsteady opposed flow flame spread over thermally 
thick solid fuels. 2 ^ All of the features were retained, 
with the exception that a simplified treatment of the gas- 
phase radiation was made, since computation otherwise would 
be too slow. 

Another current research effort 24 has examined the 
^■9 n ition and spread of a flame in a zero gravity environ- 
ment. This unsteady model is mainly different from others 
in the assumption that the velocity field can be calculated 
using potential flow, relaxing the no-slip boundary 
condition on the fuel surface. Initial results considered 
the axisymmetric case of quiescent spread. The case of 
imposing a slow flow on the flame is under development. 
Here, a three-dimensional computation is required. 

1.5. Combustion Experiments in Space 

Finally, some mention of the ultimate microgravity 
environment will be made. Ideally, we would like to have 
unlimited time, space, and accessibility to carry out 
microgravity combustion experiments. The best we can 
achieve today, with respect to duration and quality of 
microgravity, is provided by spacecraft such as the Space 
Shuttle. However because of the large amount of time and 
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money required to build and perform a space experiment, the 
study of combustion in a spacecraft environment has been 
severely limited. In fact, up until a few years ago, only 
one set of experiments had been performed. In the 1970's, 
the combustion characteristics of various practical solid 
materials were examined aboard Skylab. 25 The objectives 
were very simple, in that only visual observations were 
made on whether and how the various materials burned, and 
in some cases, extinguished, in quiescent environments. 
Motion picture photography enabled measurement of spread 
rate. In general, the flames were reported as being weaker 
than their normal gravity counterparts. However, while 
venting the chamber to extinguish the flame, combustion 
would first intensify due to the generation of air cur- 
rents. This hinted at the importance of flow on combustion 
in microgravity. 

Within the last two years, the second set of experi- 
ments, examining the combustion characteristics of thin 
paper samples, was performed aboard the Space Shuttle. 23 
Several successful runs at different oxygen percentages and 
pressures have been carried out in a quiescent chamber. 
Thermocouple readouts as well as film photography were 
utilized. Results demonstrated that the flames were in 
general weaker, more diffuse, less yellow, and larger than 
their counterparts in normal gravity. The appearance of 
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the flames suggests an extinction mechanism caused by a 
drop in flame temperature due to heat loss. This heat loss 
mechanism has been proposed in an earlier theoretical 
model . 4 

In mid-1992, several small scale combustion experi- 
ments were performed aboard the shuttle within a glovebox 
module. 26 Three experiments examining candle flames, wire 
insulation flammability, and smoldering combustion were 
successfully performed. Motion picture photography, still 
photography, and thermocouples were used. Results are 
preliminary, but the microgravity environment has once 
again produced new and interesting phenomena. 

The preceding paragraphs represent the entire history 
of spacecraft-based combustion experiments. While addi- 
tional combustion experiments are slated for spacecraft in 
this decade, by and large, most of the effort has relied on 
earth-based facilities such as drop towers, Keplerian 
trajectory aircraft, and computational studies. 

1.6. Summary 

The brief presentation above shows the current 
direction of flame spread modelling. Clearly, work 
continues to be guided by increased computational power. 
Earliest work focused on obtaining simple heat transfer 
based expressions. Then, similarity solutions were em- 
ployed. Eventually, parabolic equations were solved 
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numerically to capture the downstream region of flames. As 
computational power increased, elliptic computations first 
started to appear. The first of these assumed a flow 
field. Ultimately even the fully elliptic Navier-Stokes 
equations were solved to capture small flames or flame 
stabilization regions. 

The scope of this work then becomes evident as an 
attempt to further our understanding by studying concurrent 
flow flames with the complete equations. Elliptic equa- 
tions are used in the flame stabilization region. Down- 
stream, parabolic equations are used to capture the 
relatively long preheat region, where the hot gases heat 
the unburnt fuel. The entire preheat history of the solid 
needs to be considered, since the rate at which the fuel is 
heated affects the flame spread rate. 



Chapter 2 . Theoretical Formulation 

The problem considered in this work can be described 
as steady, concurrent flow flame spread over a thin solid 
fuel in low-speed forced flow. More specifically, "concur- 
rent" means that the flame spreads in the same direction as 
the flow velocity. "Thin" means that the fuel is both 
hydrodynamically and thermally thin. Hydrodynamically thin 
means that the thickness of the fuel is always much smaller 
than the distance of the flame from the fuel surface (or, 
equivalently, the gas phase length scale) . Thermally thin 
means that the rate of conduction in depth of the solid 
fuel is much faster than the heat up rate of the gas phase 
along the fuel surface (effectively, this means that the 
temperature is constant across the fuel thickness) . 
Furthermore, the fuel is considered thin enough that it 
burns out before the flame becomes too large (e.g. turbu- 
lent) . Low-speed forced flows are considered, at strictly 
zero gravity. 

The geometry is shown schematically in fig. 2. Here, 
the flame is shown stabilized over the fuel, which burns 
out completely at x=0 . A steady* formulation is allowed 


*For flame spread over thick fuels, the ignition 
method can affect subsequent flame spread and extinction. 
Thus, an unsteady formulation is needed. For thin fuels, 
the ignition method shouldn't matter as long as a steady 
solution exists. 
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by solving the equations in "flame-fixed" coordinates. The 
coordinate system is attached to the burnout point of the 
solid fuel. A transformation between laboratory and flame- 
fixed coordinates is given simply by considering the rate 
of flame spread. In flame-fixed coordinates, the forced 
convective stream feeds in from the left at the rate 
diminished by the flame spread rate v f, which is an unknown 
eigenvalue. In addition the fuel feeds in from the right 
at the flame spread rate. 

The model employs two-dimensional, steady conservation 
equations. The domain is divided into an "elliptic region" 
and a "parabolic region" as shown in fig. 2. These regions 
get their name from the nature of the equations being 
solved there. In the elliptic region, diffusion in both 
coordinate directions is considered and boundary conditions 
are required around the entire boundary of the region. In 
the parabolic region, diffusion in the stream-wise direc- 
tion is assumed to be small compared to diffusion in the 
cross— stream direction. Thus, boundary conditions are only 
needed on two sides of the domain (the top and bottom) . 
The parabolic region also requires an initial profile 
specification . 

One might ask why the parabolic region is necessary at 
all if the elliptic region captures the flame stabilization 
zone. The answer lies in the solution of the solid phase 
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equations . The solid begins to be heated far downstream of 
the flame stabilization zone. Therefore, it is wise to 
treat this long region with the more economical parabolic 
formulation as soon as the Reynolds number (or more 
appropriately, the Peclet number) is large enough. 

The elliptic region is considerably more difficult to 
solve both because the equations are more complicated and 
since boundary conditions are required around the whole 
region. Our goal has been to make use of the parabolic 
equations where appropriate to save computation time, while 
at the same time using the full elliptic equations to 
accurately model the flame stabilization region where 
diffusion in both coordinate directions is important. The 
two regions are described in detail below. 

2.1. Elliptic Region 

2.1.1. Flow Field Equations 

The flow field is solved using the full Navier-Stokes 
equations. The continuity equation is needed for closure. 
Stokes' hypothesis is assumed, and is a good approximation 
for air . 27 In two dimensions, the equations are: 


Continuity : 
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Conservation of momentum in y-direction: 

*f)]} 

♦WiHI)] 


(3) 


Please note that a bar above a variable indicates that it 
is a dimensional quantity. 

The inertia terms are on the left hand side of the 
momentum equations and represent the acceleration of a 
fluid element as a result of the applied forces, which 
appear on the right side of the equations. The first term 
is the pressure gradient, the second set of terms represent 
the viscous shear stress, and the last term in the x— 
momentum equation is the force due to buoyancy. This 
buoyant term takes the convenient form as shown after the 
hydrostatic pressure distribution (which would occur in a 
quiescent fluid) is subtracted. In this form, it is 
evident that a density difference leads to flow. The 
hydrostatic pressure field due to the body force term is 
not of interest. (However, for all calculations in this 
work, it is assumed that g = 0.) 

the density of the fluid changes substantially 
in a flame, the equations are written in compressible form. 
The density is evaluated by using the equation of state: 

P=pRT (4) 

The temperature dependence of viscosity is assumed to be 
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given simply as: 

jlocT (5) 

In reality, viscosity (and other transport coefficients 
presented later) are a weaker function of temperature than 
the assumed linear relation given in eq. (5) . However for 
simplicity, we use this expression, and do not expect any 
qualitative trends in our results to be affected. 

2.1.2. Flow Field Boundary Conditions 

As described earlier, boundary conditions need to be 
specified over the entire region. Please refer to fig. 2 
for the complete picture. 

Inflow ( x - x min ): u=U oc -V F ,v = 0 (6,7) 

Top (Y=Ymax): u=U„-V F , av/3y = 0 (8, 9) 

Outflow ( x = X max ) : 5u/3x and 3v/dx are given by the 

parabolic solution, which is described later. 

Bottom (y = 0) : 

x < 0 : 3"u/dy = 0 , v = 0 (10, 11) 

x > 0 : u = -V F , v = v w (12, 13) 

The spread rate is an eigenvalue which is found as part 

of the solution. It is updated at each iteration, and is 
found in the solution of the solid phase equations, 
described in a later section. The velocity v w is also 
determined by the solid phase equations. Basically, the 
mass flux from the fuel depends only on the temperature of 
the fuel. Knowing the mass flux, we can calculate v w . 
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Since the velocity is given at the inflow and bottom 
of the domain, pressure is not known here. Because very 
low speed flows are considered, pressure varies only 
slightly from the ambient value of 1 atm. For this 
problem, it is sufficient to set pressure to 1 atm. at an 
arbitrary point in the domain. The upper left corner is 
chosen . 

These boundary conditions assume that the flow is 
perfectly aligned with the fuel plate. While in a space- 
craft environment the recirculation currents are quite 
random, the sensitivity of the flame to angle of flow 
velocity was outside the scope of this effort. 

Through these boundary conditions, one of the princi- 
pal parameters of interest is varied, namely the free 
stream velocity U,- . The last boundary conditions, eqs . 
(12) and (13), are given by the solution of the solid phase 
equations. From the point of view of the gas-phase 
equations, the coupling between the gas and solid phases is 
made through these boundary conditions. 

2.1.3. Nondimensional Equations and Boundary Conditions 

There are several reasons for using nondimensional 
equations. A proper choice of scaling parameters simpli- 
fies numerical computation. Consider an example where we 
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use a length scale*, L 0 , to nondimensionalize the spatial 
coordinates in the governing equations. We want to ensure 
that L 0 is the smallest scale of importance in the physical 
solution of the problem, so that we know that taking ten 
grid points in that scale, for example, should be suffi- 
cient to accurately solve the problem. Thus, much of the 
guesswork in choosing a grid is eliminated as other 
parameters in the problem are changed. Another important 
reason to nondimensionalize is to capture the important 
parameters in the problem, for example, the Reynolds 
number. Finally, from a practical point of view, nondim- 
ensional variables are good to work with since their values 
should be order unity. Otherwise, it is sometimes clumsy 
to deal with very small or very large numbers in the same 
problem, and can make the presentation of results in such 
cases somewhat confusing. 

Velocity is nondimensionalized by the characteristic 
relative velocity Ur = U,„ + U B - V F . The three components 
represent a contribution from forced flow, buoyant flow, 
and flow due to the rate of flame spread, respectively. u r 
is a measure of the net velocity near the flame stabili- 
zation region due to these three terms. 


*The length scale is a function of other parameters 
in the problem, and thus can change with these parameters 
for different cases. 



30 


The coordinates ic and y are nondimensionalized by the 
thermal length, given as x R = a */U R . The thermal length is 
found by considering a conduction-convection balance, and 
is a good measure of the flame standoff distance and the 
flame thickness in the stabilization zone. The asterisk 
indicates that the property is evaluated at the temperature 
T * , which is the mean of the adiabatic flame temperature 
and the ambient temperature. For convenience, we use 
T* = 1250 K in all cases. 

The property values are based on air, which is a good 
^PP^oxiroat ion to the dominant component in the gas phase. 
While certain property values depend greatly on composition 
(e.g., large heavy fuel molecules behave quite differently 
than air) , a detailed specification of property values 
based on composition is not attempted at this time. One 
reason for this is that the chemistry assumed is rather 
simplified, so any prediction based on the chemistry would 
be very limited. 

Density and viscosity are nondimensionalized by their 
values at the reference temperature T * . Pressure is 
referenced to the ambient value of 1 atm. and is nondim- 
ensionalized as p = (P-Poo) /P*U r . These (and all) non- 
dimensional parameters are also listed in the nomenclature 
list. Using these quantities, the equations become: 
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Continuity : 


d(pu ) + d(pv ) _ 0 
3(x) 3(y) 


(14) 


( 15 ) 


(16) 


Conservation of momentum in x-direction: 

* Pi £K§7*£)] tGr(p *' p> 

Conservation of momentum in y-direction: 

In these equations, the Prandtl number, Pr, appears. This 
is the result of choosing the thermal length as the 
characteristic length scale. For all cases, Pr is assumed 
constant and equal to 0.7. The Reynolds number based on x, 
Re x , is given simply as x/Pr, and similarly, Re y =y/Pr. 
Equivalently, x and y represent Peclet numbers based on x 
and y respectively. The boundary conditions become: 

Inflow (x = x min ) : »- (U.-V F )/U„, v = 0 (17,18) 

Top (y = y max ) : u = (U~ - / ^ v /^y = 0 (19/ 20) 

Outflow (x = x max ) : du/3x and 9v/8x are given by the 

parabolic solution, which is described later. 

Bottom (y = 0) : 


x < 0: 3u/dy = 0 , v = 0 (21, 22) 

x > 0: u = -V F /U R , v = v w (23, 24) 

Note that for purely forced flow, (U~ - V F ) /U R = 1 . 

2.1.4. Energy Equation 

By assuming the specific heat of the gas is constant. 
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conservation of energy enables us to write an equation for 
temperature as follows: 


PuC>-g*P^.f 



(25) 


The left hand side of the equation represents the convec- 
tion of heat. The right hand side is the heat conduction 
term and the source term due to the chemical reaction. The 
heat release due to viscous dissipation and compressive 
work are neglected since they are small compared to the' 
combustion heat release. The specific heat, c p , is assumed 
constant. The conductivity is assumed to vary according 
to : 


K«T (26) 

The form of the energy source term is: 


Q =Q B fl p 2 Y F Y 0 exp(-E ff /RT) (27) 


This finite rate expression is for a one-step, second order 
reaction . 

2.1.5. Energy Equation Boundary Conditions 

The boundary conditions for the energy equation 
resemble those for the momentum equation, since they are 
both elliptic in nature: 

Inflow (x = x min ) : T^ (28) 

Top (7 = ymax): T = T~ (29) 

Outflow ( x = X max ) : 9t/8x is given by the parabolic 

solution, which is described later. 



33 


Bottom ( y =0 ) : 


x < 0 : 3T/3y = 0 (30) 

x > 0 : T = T s (31) 


The relation given in eq. (31) is determined by solution of 

the solid phase equations, which is described shortly. 

2.1.6. Nondimensional Energy Equation 
and Boundary Conditions 

Temperature is nondimensionalized by the ambient 
temperature, L, and conductivity is nondimensionalized by 
its value at T* (of these two reference temperatures, the 
former is useful in presenting results, while the latter is 
appropriate for evaluating property values) . The equation 
becomes : 


PU ^ 


d r 

tpv 0y 



(32) 


where 

Q = DaQp 2 Y F Y 0 exp(-E/T) (33) 

Note the appearance of the Damkohler number, Da, in the 
last equation ( Da = a* p* B g /Ur ) . The Damkohler number is 
the ratio of the characteristic flow time over the chemical 
reaction time in one thermal length. It varies inversely 
as the square of relative velocity. The boundary condi- 
tions are: 

Inflow (x = x min ) : T = 1 (34) 

Top (y = y max ) : T = 1 (35) 

Outflow (x = x max ) : dT/3x is given by the parabolic 
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solution, which is described later. 

Bottom (y = 0) : 

x < 0 : dT/dy = 0 (36) 

x > 0: T = T s (37) 

In eq. (37) , T s is determined by solution of the solid 
phase equations. 

2.1.7. Fuel and Oxygen Species Conservation Equations 

The species conservation equations are: 
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The left hand side of the equations represents the convec- 
tive term. The right hand side is the mass diffusion term 
and the sink term due to the chemical reaction. The 
diffusion coefficient, P D i, is assumed to vary as follows: 

(pD^ocT (40) 


The form of the species sink terms is: 

= -B g p 2 Y F Y 0 exp(-E g /RT) , <5^ = f<5; (41, 42) 


where f is the stoichiometric oxidizer/fuel mass ratio. 
2.1.8. Fuel and Oxygen Species Boundary Conditions 

Again, the boundary conditions resemble those in 
preceding sections since the equations are elliptic: 
Inflow (x = x min ) : Y F = 0, Y 0 = Y 0<oo (43, 44) 



35 


Top ( Y = Ymax ) : Y f = 0, Y 0 = Y 0>oo (45, 46) 

Outflow ( x = X max ) : 0Y f /0x and 0Y o /dx are given by the 

parabolic solution, which is described later. 

Bottom ( y =0 ) : 

x<0: 0Y F /0y = O, 0Y o /0y = O (47, 48) 


mY F =in + pD F (0Y F /3y) w (49) 

x>0: ffiY 0 ; v = pivavay). <so> 

2.1.9. Nondimensional Species Equations 
and Boundary Conditions 

The terms P D i are nondimensionalized by their value at 
T*. The equations become: 




Where: cb F = -Dap 2 Y F Y c exp (-E g /T) , tb 0 = fci> F (53, 54) 


The boundary conditions are: 


Inflow (x = x min ) : 

Y F = 0, Y 0 = Y 0> „ 

(55, 

56) 

Top (y = y max ) : 

Y f = 0, Y 0 = Y 0fOO 

(57, 

58) 


Outflow (x = x max ) : 0Y f /0x anc i 0Y o /0x are given by the 

parabolic solution, which is described later. 

Bottom (y = 0) : 

x < 0: 0Y F /0y = O, 0Y o /0y = 0 (59, 60) 

^ > rhY F , w =m + (pD F /Le F ) (0Y F /0y ) w (61) 

X ” mY 0(W = (pD 0 /Le 0 ) (0Y o /0y) w (62) 
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2.2. Parabolic Region 

The equations for the parabolic region are identical 
to those in the elliptic region except for the neglection 
of the stream— wise diffusion terms. In addition, in 
deriving the flow field equations, the y-momentum equation 
drops out completely, and, in this case, the pressure 
gradient term in the x-momentum equation is zero. The 
equations have the well-known form of boundary layer flow 
over a flat plate. They are presented below. 

2.2.1. Nondimens ional Parabolic Region Equations 

Continuity: ^TT + = ° <63) 

X-Momentum: pu-^ + p v -|^ = Pr ) +Gr “P^ (®4) 

Energy: P u j£ + P v (65) 
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2.2.2. Nondimensional Parabolic Region Boundary Conditions 

The parabolic region requires boundary conditions only 
on the top and bottom of the domain, as well as an initial 
profile specification. The initial profile is given by the 
values obtained at the outflow (x = x max ) of the elliptic 
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solution. Thus, velocity, temperature, and the species 
concentrations are given. The boundary conditions on the 
top and bottom of the domain are identical to those 
presented in previous sections.* 

2.3. Coupling Between Elliptic and Parabolic Regions 

The coupling between the elliptic and parabolic 
regions is best understood by considering the solution 
algorithm. First, the elliptic equations are used for 
several iterations to get an updated elliptic field. Then, 
the values of the variables (u, v, T, Y F , and Y 0 ) at the 
outflow of the elliptic region are used as the initial 
profile for the parabolic equations. After sweeping 
through the parabolic region and the solid phase equations, 
the elliptic region is again calculated. This time, 
however, the boundary conditions at the outflow for the 
elliptic equations are based on the parabolic solution, by 
evaluating derivatives with respect to x in the parabolic 
region at the interface between the two regions, and using 
these derivative boundary conditions for the elliptic 
region . 

In previous elliptic computations of this sort, the 
outflow boundary conditions are given simply as a zero- 

*One exception is the boundary condition for v-ve- 
locity. Since the y-momentum equation drops out, v needs 
to be specified only on the bottom of the domain. 
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gradient normal to the boundary (e.g., ref. 14) . This is 
because information beyond the boundary is unknown. A 
result is that the solution near the boundary is not as 
accurate as the solution away from the boundary. In our 
case, we can use the actual (non-zero) gradient as given by 
the first step of the parabolic computation, as described 
above. Upon iteration, the converged solution provides a 
much smoother transition between the elliptic and parabolic 
domains . 

This is shown in fig. 3. In fig. 3a, zero gradient 
boundary conditions for the outflow of the elliptic region 
are used. In the elliptic region, the isotherm near the 
boundary levels off due to the zero-gradient condition. 
Thus, a slight dip is evident. 

In fig. 3b, the results of the method used in this 
work are shown. Now, there is a smooth transition between 
the two regions. 

2.4. Solid Phase Equations 

The solid phase is solved by considering the conserva- 
tion of mass and energy together with a pyrolysis relation. 
The flame spread rate, v f, appears as an eigenvalue. The 
geometry of the solid fuel is shown in fig. 4. The density 
of the solid is assumed constant while the thickness is 
allowed to decrease due to pyrolysis. In fig. 4 dimensions 
are greatly exaggerated. As far as the gas phase is 
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Figure 3. Elliptic / pa ra bolic matching boundary conditions. 
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Figure 4. Geometry of Solid Fuel. 
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concerned, the solid is infinitesimally thin, meaning that 
the flame stand-off distance is much larger than the fuel 
thickness. In deriving these equations, it is assumed 

that : 

y 1 + (dh/dx ) 2 ~1 (68) 


This is an excellent assumption for all cases studied 
The solid phase equation is: 


— — — d(hT s ) _ /rp4 t M 

Qy + P 3 V F C s ^ °® < T S 2 Tm ' 


+ P S V F -^C( C B - C p) T L + (- L) +C P T sl (69) 


On the left side, the first term, <* y , is the conductive 
heat flux incident on the solid from the gas phase. It is 
given as q y = t K ( 3 t/ 3y ) ] g , w . This term represents the 
coupling between the gas and solid phase from the point of 
view of the solid equations. The solid phase spans both 

the elliptic and parabolic regions. The second term 

represents the convection of heat due to flux of fuel. As 
described earlier, in the flame-fixed coordinate system 
used, the fuel feeds into the domain at the speed v f . It 
may be clearer for the reader to examine this term in 
laboratory coordinates. In laboratory coordinates, it is 


simply the unsteady bulk heat up term. 

On the right side, the first term is the surface 
radiative heat loss term. The second term is the energy 
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change due to the latent heat of vaporization of the fuel. 

L is the latent heat of the fuel, which is specified at the 
temperature T L (T l = 300 K in this problem) . 

The pyrolyzed gases constitute blowing from the wall. 
The pyrolysis law chosen specifies that the rate of pyroly- 
sis depends only on the temperature of the fuel (zeroth- 
order pyrolysis) . it is: 


ffi= A e p 3 exp(-E a /RT 3 ) (70) 


Using eg. (70), the blowing velocity v w can be foundf since 
the density of the gas is known at the wall. Remember that 
the solid density, P s , is constant. This is different from 


other models 11, 15,22 


which assumed the thickness of the fuel 


was constant, but the density was decreased as pyrolysis 

occurred, making the pyrolysis relation first-order. The 

reason for using a zeroth-order pyrolysis relation is that 

complete fuel burnout is possible. A first-order pyrolysis 

law requires an arbitrary (but non-zero) specification of 

unburnt solid char, at which point burnout is said to 
occur . 


When combined with the conservation of mass for the 
solid fuel, the pyrolysis relation becomes: 

dh A — — 

— = exp(-E 8 /RT 8 ) (71) 

V p 

Equations (70) and (71) are solved with the following 
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boundary conditions: 

At x = x maX(par : T=T„, h=T (72, 73) 

The flame spread rate is an eigenvalue and is solved for 
iteratively. In solving the solid phase equations, the 
conductive heat flux distribution from the gas phase is 
known as a function of x. The spread rate, 7 F( is guessed. 
If this guess is too high, the fuel moves too fast for 
complete burnout to occur. If the guess is too low, 
burnout occurs too soon (i.e. at x" > 0 ) . The spread rate is 
adjusted using a simple bisection method until burnout 
occurs at x = 0 . 

Using x R/ p'*CpU R T 00/ and C s to nondimensionalize h , 
q y , and C p respectively, the equations and boundary 
conditions become: 


+ E d(hT) = r ( T ^ -1) (1 -C p ) T L + <-L) + Cp T ] 


At x = x 


max, par * 


T = 1, h = t/x f 


(75, 76) 


In eq . (74), >— - Ps^f^s ! P* / and r-0£T w /p /C P /U R . 

The last term, T, is the non-dimensional surface radiative 
heat loss parameter. It varies inversely with relative 


velocity . 



Chapter 3 . Numerical Model 

The problem is analyzed by writing the finite-differ- 
ence representation of the equations and solving them using 
the CRAY X-MP supercomputer at NASA Lewis Research Center. 
The algorithm used is called SIMPLER (Semi-Implicit Method 
for Pressure-Linked Equations Revised) . This algorithm is 
presented in detail in a textbook, and so only the unique 
features are presented here. The computer program, which 
was written based on this algorithm specifically for this 
problem, appears in the appendix, along with parameter 
values . 

In any finite-difference method, the domain is broken 
up into a grid structure. The equations are written for an 
individual control volume (or grid square) by considering 
its interaction with neighboring control volumes . The 
solution can be thought of as equivalent to finding the 
parameters of interest at the interfaces of adjacent 
control volumes by an appropriate interpolation method, and 
then calculating the quantities at all control volumes 
using some matrix inversion technique. Finite-difference 
methods vary according to how they determine interface 
quantities. Simple interpolation methods require less 
computational complexity, but suffer in accuracy. 

As can be seen in the Theoretical Formulation chapter, 
the equations are quite complex. Both first and second 


44 



45 


derivative terms appear. In addition, the equations are 
non-linear. These attributes lead to several complica- 
tions . 

3.1. Discretization of Convective Terms 

The terms such as pu(8T/3x) representing convection, 
are somewhat troublesome. The problem lies in writing the 
appropriate finite difference representation for a given 
control volume. The physical meaning of such a term should 
be considered when deriving the finite difference form. If 
we assume that temperature, T, represents the energy of the 
fluid, then pu(3T/3x) is the rate of energy transfer to 
the control volume due to the fluid flowing at velocity u. 
This is of course why it is called a convective term — the 
energy is convected by the velocity. Clearly., anything in 
a stream of fluid is generally affected more by the 
upstream than the downstream flow. In writing the discret- 
ization equations, then, it is important to weight the 
upstream quantity. It is the directional character of 
velocity which drives this consideration. One common 
approach is to use the upwind-difference scheme (or donor- 
cell method) . In this method, the value of a variable 
(say, temperature) at the interface of a control volume is 
set equal to the value on the upwind side of the inter- 
face. This is perhaps the simplest possible weighting 


function . 
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While the upwind-difference scheme captures the 
essential elements of a convective-diffusive transport 
problem, there are instances when the method becomes less 
accurate. Specifically, as the grid Peclet number*, Pe, 
becomes large, conduction is overestimated. This is 
because while at an interface, the value of temperature is 
correctly given as described above, the conduction is 
assumed based on a linear profile between the two points on 
either side of the interface. In actuality, the conduction 
at the interface is close to zero, but the upwind scheme 
predicts a value much larger. 

To overcome this difficulty, the exact solution of the 
convective-diffusive equation 



is determined for adjacent control volumes. This solution 
determines the value of temperature, T, to use at the 
interface of the control volumes. Consider two adjacent 
control volumes separated by Ax . Let T L and T R be the 
temperature of the two control volumes, the exact solution 
of eq. (77) (assuming pu is constant for the control 
volumes) is then: 


*The grid Peclet number is a measure of the relative 
strength of convective to conductive transport in an 
individual control volume. 


( T (x) -T l ) _ [exp(puC p x/ K) - l] 
(Tr -T l ) texp(puC p Ax/K) - 1] 


( 78 ) 


and the value at the interface between the two control 
volumes is readily given. Since exponential terms are 
relatively costly to compute, an approximating polynomial 
(fifth order in x) is used instead. The approximation is 
very close to the exact solution, enabling accuracy for any 
grid Peclet number. 

While eq. (77) is one-dimensional, the same idea is 
used to solve the two-dimensional equations in this work. 
The convective-diffusive transfer is simply considered 
along each coordinate direction. 

3.2. Treatment of Source Terms 

In the SIMPLER algorithm, all terms other than 
convective and diffusive ones are lumped together and 
called source terms* . The rules governing their discret- 
ization follow. 

Suppose the equation of interest is the conservation 
of energy, for which temperature is the dependent variable. 
The source term is in general a function of T, derivatives 
of T, and other variables in the problem. The key is to 


*The only exception to this is in the solution of 
the Navier-Stokes equations where the pressure gradient 
terms are handled individually. This will be described 
subsequently . 
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linearize the source term as 

S = S c + S p T (79) 

and ensure that S p is always less than zero. In eq. (79), 
T represents the temperature of a given control volume, 
which is to be calculated. The coefficients S c and S p may 
be functions of any number of variables in the problem 
(including a value of T that is guessed or taken from a 
previous iteration) . This specification guarantees that 
the source term will not lead to instability, even for the 
highly non-linear Arrhenius rate expressions of combustion 
in this problem. 

3.3. Staggered Grid System 

In developing the finite difference form of many 
equations, it is generally a good idea to use the same grid 
structure for each dependent variable to minimize the 
complexity in bookkeeping all of the grid locations and 
interpolation between grids. However, there is nothing 
wrong with using different grid structures for different 
equations. It only makes sense to do this if some benefit 
can be derived. In the SIMPLER algorithm, a staggered grid 
is used to eliminate problems brought on by the appearance 
of first derivative terms, such as pressure gradient terms 
or terms appearing in the continuity equation. 

A very coarse grid structure is shown for demonstra- 
tion purposes in fig. 5. The grid is generated by first 


Main Variables U-velocity 
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and u and v — velocities. 
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chopping up the overall domain into smaller control 
volumes. The grid location for each dependent variable 
(including pressure) is taken to be the center of the 
control volume. Then, u-velocity is positioned at each 
vertical face, and v-velocity at each horizontal control 
volume face. Thus the composite grid is comprised of three 
different grid structures. 

The advantage of doing this is that in evaluating the 
pressure force on a velocity control volume, the pressure 
at adjacent grid points is used. If all variables (includ- 
ing velocity) were evaluated at the same grid points, the 
pressure force would be calculated by using every other 
grid point . This would be the result of a simple discret- 
ization of the first derivative pressure gradient terms. 
Unless treated specially, this kind of a formulation 
permits highly unrealistic solutions. 

3.4. The SIMPLER Algorithm 

To obtain a solution, the initial variable field is 
taken from a converged case. Then, the velocity or oxygen 
percentage is changed by modifying the boundary conditions. 

The first step of the algorithm is to start with the 
non— converged velocity field, which is supplied initially, 
or taken from the previous iteration. Then, the momentum 
equations are used to obtain an expression for the velocity 
(x and y-direction) of the fluid in a control volume. 
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These finite difference expressions will of course include 
pressure. By substituting the expressions for velocity 
into the finite difference form of the continuity equation, 
an equation for pressure results. Pressure (call it P*) is 
calculated*. Using this updated pressure field, the 
momentum equations are solved to obtain new estimates of 
velocity . 

During iteration, these velocity estimates will not 
satisfy the continuity equation. To speed convergence, a 
superposed pressure field, called P', is assumed. Thus, 
P=P*+P', and after convergence is reached, P'=0 everywhere. 
Using the momentum equations, simplified relations for 
velocity involving P' are obtained. These are substituted 
into the continuity equation to obtain a P' equation. 
After computing P' , it is used to correct the velocity 
components . 

At this point, all additional conservation equations 
are solved (in this case, conservation of energy and 
species) . Finally, the entire procedure is repeated until 
convergence is reached. 


‘in the SIMPLER algorithm, the equations are dis- 
cretized such that a tri-diagonal matrix results for a 
given line of grid points. The tri-diagonal matrix is 
easily inverted. Thus, in two dimensions, each row of 
grid points is solved by matrix inversion, followed by 
each column of grid points. This procedure pulls the 
boundary information into the domain very quickly. 
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3.5. Domain Structure 

The grid structure should be fine enough to adequately 
capture the changes in the parameters of interest. For 
example, in a flame spread problem, temperature varies 
dramatically across the flame. Thus many grid points are 
required. However, if a uniform distribution of grid 
points is used throughout the domain, most of the grid 
points are unnecessary in regions where the temperature 
does not change as dramatically as in the flame. Therefore 
it is more efficient to use a variable grid spacing, i.e., 
in regions where the quantity varies sharply, a higher 
concentration of grid points is used.* 

For the flame stabilized at the leading edge of a fuel 
plate, the grid points should be concentrated at the 
leading edge both to capture the high temperature, narrow 
flame stabilization region as well as the rapid flow field 
changes due to the sudden appearance of the plate. The 
variable grid structure used in this problem is shown in 
fig. 6. Evident in fig. 6 is the concentration of points 
both immediately upstream and downstream of the leading 
edge to capture the flame stabilization zone and the abrupt 
pressure and velocity changes. Additionally, the concen- 


*An adaptive grid method or multigrid method would 
be well-suited for this problem, but are not used at this 
time . 
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tration of points near the fuel plate capture the steep 
gradients near the plate. The location of a typical flame 
is shown superimposed on the grid structure. The particu- 
lar grid spacing used in this problem was found to be more 
than adequate in an earlier work. 14 

In the elliptic region, x varies from about -20 to 70 
and y varies from 0 to 90. The number of grid points in 
the x-direction is 70 and the number of grid points in the 
y-direction is 50. 

The Reynolds number*, Re x , as a function of x, for any 
case is simply calculated as x/Pr, where Pr is the Prandtl 
number (Pr=0.7). For all cases in this effort, the 
elliptic/ parabolic boundary occurred at Re x =100. A test 
case was run where the elliptic/parabolic boundary was 
extended to Re x =200. The thermal structure and spread rate 
of the flame changed by less than 2%, so the smaller 
elliptic domain was deemed adequate. 

3.6. Parabolic Region 

The finite difference scheme in the parabolic region 
is based on a marching technique. 29 At the outflow of the 
elliptic domain (Re x =100 or x=70), parabolic equations take 
over. The results obtained from the elliptic region 
provide the initial condition for the parabolic equations. 

*Re x is based on the relative velocity, distance from 
the leading edge, and properties evaluated at T*. 
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A single sweep into the far downstream region (Re x =2200 or 
x=1570) is utilized. 

The same, variable grid structure in the y-direction 
as used in the elliptic region is utilized. Also, addi- 
tional points are added to the top of the domain (at y>y max 
in fig. 2) to give the boundary layer room to grow. Thus, 
in the parabolic region, x varies from 70 to 1570 and y 
varies from 0 to 150. The grid spacing in the x-direction 
gradually increases, since derivatives with respect to x in 
the far downstream region become smaller and smaller, 
making a coarser grid adequate. There are 65 grid points 
in the x-direction and 90 grid points in the y-direction. 

An implicit formulation is employed, that is, at any 
station x=x ir the resulting equations depend only on the 
values of the variables at x=x i _ 1 (which are known) and at 
x=x i . To get the solution at x=x ir a tri-diagonal system of 
equations is solved. 

3.7. Computer Usage 

To compute a new case, a nearby, converged flame. is 
used as the initial guess. One global iteration step takes 
about 7 sec. of CRAY X-MP CPU time. It consists of first 
calculating the elliptic region using line-by-line sweeps 
(4 horizontal sweeps, 4 vertical sweeps) for each of the 
conservation equations in series. This is repeated 10 
times. Then, the parabolic region is calculated using the 
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marching technique described above. Finally, the conduc- 
tion of the gas phase to the solid is used to solve the 
solid phase equations and get new values of spread rate, 
temperature, and blowing velocity. Using these values, the 
boundary conditions are updated, and the global iteration 
step is concluded. Typically, convergence requires ab out 
1 hr. of CPU time. 

3.8. Property Values 

The property values used are presented in the nomen- 
clature list. They are based on the choice of cellulose as 
fuel. Some of the derived parameters require description. 

Consider cellulose and air stoichiometric combustion: 
C 6 h ic°5 + 6 (0 2 + 3.76N 2 ) 6C0 2 + 5H 2 0 +22.56N 2 (80) 

Assuming one mole reacts: 

( 1 62 g) c 6Hlo05 + (192 g) 02 + (6 3 2 g)^-* 

(264 g) COz + (90 g) h 2 o + (632 g) N2 (81) 

In order to calculate p*, the mixture specific ideal gas 

constant used is the average value of that for air and for 

the products (nitrogen is included in the products) : 

R = ( R Air + R Procj) / 2 = (0. 287 + 0. 302 )/2 = 0. 295 J/g/K (82) 

Based on the assumed adiabatic flame temperature of 
2200 K, the gas-phase specific heat is calculated from the 
simple energy balance as follows: 

(Q + L ) m Fuel = m ProdCp (T f - T^) (83) 

and since iti Prod /m Fuel = 6 . 084 f W e get C p = 0.30 cal/g/K. 



Chapter 4 . Results 


Concurrent flow flame spread over a thin fuel in a 
zero gravity* environment has been modelled. Steady 
computations have been carried out over a range of forced 
flow velocities and ambient oxygen concentrations at one 
atmosphere pressure. The molar* oxygen percentage was 
varied from 13.5 to 21, and the flow velocity was varied 
from 16 cm/s down to 0.8 cm/s*. As will be described 
later, the lower range of both parameters was determined by 
an extinction boundary. The upper ranges were chosen arbi- 
trarily, but were high enough to adequately present the 
desired trends. The property values, which are given in 
the nomenclature list, correspond to a thin cellulosic 
sheet, which has been used extensively in recent micro- 
gravity experiments? 0,30 ' 31 

While the computation was carried out in nondimension- 
al coordinates, some results will be presented in dimen- 
sional form. As a reminder, the space coordinates (x and 
Y) were non-dimensionalized by the thermal length. The 

*Low-speed forced flow at strictly zero gravity is 
modelled. The effect of microgravity is not considered at 
this time. 

*Unless otherwise noted, the oxygen content will 
always be expressed on a molar basis. 

*The magnitude of all of the forced flow velocities 
considered were well below those encountered in normal 
gravity flame spread. 
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elliptic domain extended from x=-20 to 70 and y=0 to 90. 
The parabolic domain extended from x=70 to 1570 and y=0 to 
150. Over the range of velocities in this study, the 
thermal length varied from 2.6 to 0.15 cm. 

A total of 40 cases have been found. The computation- 
al test matrix is shown in fig. 7. Each circle on this 
plot represents a converged, steady flame. The two 
parameters varied were flow velocity and oxygen percentage. 
4.1. Detailed Flame Characteristics For One Case 

Before proceeding with any comparisons, one case will 
first be presented in detail which shows all of the 
features that the computational model is capable of. For 
this flame, the oxygen content was 15% and the free stream 
velocity 5 cm/s. The thermal length for this case was 0.46 
cm. This typical case was roughly in the center of the 
parameter space. 

4.1.1 Gas Phase Profiles; 5 cm/s, 15% 0 2 

Fig. 8 shows isotherms in the nondimensional space. 
The temperature field is able to be probed extensively. In 
part (a) , a long view is shown to demonstrate the extent of 
the flame and thermal plume. This view in fact covers a 
portion of the elliptic and parabolic domain. More will be 
said later, but it is evident that the isotherms smoothly 
span the interface (x=72) between the two domains. The 
maximum temperature in the gas phase occurs very near the 
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burnout location (x=0) , as indicated. This is true for 
most of the cases examined. Only as the flame approaches 
blowoff does the maximum temperature get shifted down- 
stream. In part (b) , the flame structure at the leading 
edge is enlarged to show the detail of the heat transfer in 
this region. Clearly, the diffusion is two-dimensional 
since significant heat and mass transfer upstream of the 
fuel has occurred. This is one of the justifications for 
using an elliptic treatment to capture the flame stabiliza- 
tion region. Notice that the temperature changes from the 
ambient value to the maximum value in the span of about one 
or two units. This verifies the choice of thermal length 
as the important length scale. In part (c) , the solution 
reflected about the y=0 plane to show what the thermal 
structure for the entire flame leading edge looks like. 
The temperature tends to fan out somewhat as it cools. 
This is a characteristic of the small, short flames 
predicted in this work. In contrast, a flame burning in a 
normal gravity buoyant environment tends to stay closer to 
the fuel in the downstream region. This flame shape effect 
is demonstrated in the small, low-velocity flames of refs. 
20, 30, and 31. 

In fig. 9, the mass transfer aspects of this diffusion 
flame are presented. In part (a), the fuel and oxygen mass 
fraction contours are shown in the leading edge region. 
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The contours overlap most near the burnout point. Notice 
how slowly the fuel is consumed in the downstream region, 
based on the large separation of contours there. In part 
(b) , the local equivalence ratio* contours are shown. 
Again, the contours are drawn out downstream. In some 
earlier works, the flame length is defined as equalling the 
length of the <|>=1 contour. 15 This definition is an artifact 
of the flame sheet models which defined the flame as 
existing at the location of stoichiometric equivalence 
ratio . 

In examining the expression for reaction rate, it is 
seen that the three requirements needed for a robust 
reaction are: 1) presence of fuel, 2) presence of oxygen, 
and 3) an elevated temperature. This common sense descrip- 
tion suggests that the maximum reaction rate for this case 
should also be confined to the near-burnout region since 
this is where the temperature is highest and significant 
oxygen and fuel overlap occurs. Indeed, in fig. 10, part 
(a), the reaction rate contours are clustered near x=0. 
The highly non-linear Arrhenius expression for the reaction 
rate is unmistakable in the very sharp gradients (the 
contours are separated by a factor of ten) . 


*The local equivalence ratio is defined as the ratio 
of fuel to oxygen, divided by the stoichiometric fuel to 
oxygen ratio. 
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One way of defining the visible flame is to choose a 
reactivity contour. This is perhaps the best way in that 
the reactivity level should roughly scale with the number 
of photons emitted in the reaction such that our eyes can 
register something there. In this work, the reactivity 
contour equal to 1CT 4 g/cm 3 /s is arbitrarily defined as the 
boundary of the visible flame (as seen later, this gives a 
flame shape very similar to that found in experiment) . In 
a numerical study such as this, it is important mainly to 
use a consistent definition of the visible flame. For 
comparison, part (b) gives the isotherms and 0=1 contour on 

the same scale. 

A benefit of solving the Navier-Stokes equations is 
that we can probe the velocity field in detail. In fig. 
!1, velocity vectors are plotted in the flame leading edge 
region*. There are several features which should be 
pointed out. As the flow approaches the fuel plate (x=0) , 
it decelerates and is deflected up. Initially, the flow 
decelerates due to the plate and the hot flame (the visible 
flame, as defined above, is shown) . The flow is deflected 
both merely by the presence of the plate and due to the 
influence of the blowing of pyrolysis products from the 


*The velocity vectors in 
ized by the relative velocity, 
cm/ s . 


fig. 11 were nondimensional- 
which for this case was 4.62 
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Figure 1 1 . Velocity vector plot. 14=5 cm/s, X 0 =1 5% 
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fuel. Since there is gas expansion due to the heat 
release, the flow accelerates in the flame zone. As we 
move downstream, the velocity profile flattens out due to 
viscous effects and the cooling of the gas stream. As 
described earlier, the equations have been solved in flame 
fixed coordinates. The fuel feeds in from the right at the 
flame spread rate. In fig. 11, the velocity on and near 
the fuel plate is indeed to the left. If the transforma- 
tion is made back to the laboratory reference frame, the 
velocity vectors in fig. 12 result. Here, the flow 
characteristics near the plate are quite conventional. The 
fuel issues from the surface in the y-direction and is 
blown downstream. 

In fig. 13, streamlines in both flame and laboratory 
fixed coordinate systems are shown. The stagnation 
streamline is chosen as ^O. It shows the influence of the 
fuel blowing, as it is lifted off the surface. The 
streamlines given by 'FcO are a result of pyrolysis gases 
being blown off the fuel surface. The streamlines given by 
¥>0 represent the imposed forced flow. 

In the velocity and streamline plots for this case, 
the difference between flame and laboratory fixed coordi- 
nate systems is small because the spread rate is a small 
percentage of the relative velocity. The spread rate is at 
most 18% of the relative velocity in all cases studied. 
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Figure 1 2. Velocity vector plot (lab. coordinates). l!*=5 cm/ s, Xq— 1 5% 
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The pressure field in the flame leading edge region is 
depicted in fig. 14. Remember that nondimensional pressure 
is used, referenced to ambient pressure, p= ( p-p oo) /p /U R . 
The magnitude of the pressure changes in dimensional terms 
is tiny, because of the very low speed flows considered. 
Two pressure rises, the first due to the presence of the 
hot flame and the second due to the plate are shown. 
Downstream, the pressure slowly returns to the ambient 
value as the gas flows from the influence of the leading 
edge and flame. 

4.1.2 Solid Phase; 5 can/s; 15% 0 2 

The solid phase will now be presented. In fig. 15, 
the heat flux incident on the solid, the solid temperature, 
thickness*, and blowing velocity are all shown over the 
whole domain. All variables are in nondimensional terms. 
The heat flux curve is quite spiked at the burn out point. 
This is because the flame gets closest to the fuel there. 
The surface temperature is relatively level in the pyroly- 
sis zone. At about x=8, there is an inflection point in 
both the solid temperature and heat flux curves. This 
indicates the end of the pyrolysis region. At this point, 
the solid thickness is at 99.9% of its unburned value. 
This criterion for solid thickness is used to find the 

*Here, thickness is nondimensionalized by its initial 


value . 
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pyrolysis length in all cases (used later in §4.3) , and 
corresponds very nearly to the inflection points evident in 
the T s and q w curves. Most of the fuel blowing into the gas 
phase is confined to a small region on the x-axis. 

4.2. Comparison With Experiment 

In ref. 30, a number of experiments were performed in 
the 5.18 sec. Zero Gravity Facility at NASA Lewis Research 
Center in Cleveland. The effort focused on concurrent flow 
flame spread over a thin fuel (tissue paper, Kimwipes*) in 
very slow speeds. The forced flow was generated by moving 
the entire fuel plate into a quiescent oxidizer. 

While the fuel sheet was not perfectly flat or uniform 
in appearance, care was taken to get as good a sample as 
possible. The detailed modelling of the macrostructure and 
pyrolysis of the fuel was outside the scope of this effort. 
We use a simplified model to capture the important trends. 

For most of the cases, 5.18 seconds was not enough for 
the flames to reach steady state. However, some of the 
flames seemed to be close to steady. One such case will be 
used for comparison. The experiment was performed at 15% 
0 2 and 4.84 cm/s. The theoretical comparison is done at 
15% 0 2 and 5 cm/s. The experiment recorded only the 
visible flame using motion picture photography. Therefore, 
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only a comparison of the overall flame dimensions and 
spread rate is possible. In fig. 16, the computed fuel 
reactivity contours are compared to the reported experimen- 
tal flame shape (this is how the flame looks just prior to 
the end of the 5.18 sec. drop). The flame shape is 
predicted quite well if the proper reactivity level is 
chosen to define the threshold of visible emission. For 
this case, that level is around 10 -4 g/cm 3 /s. In the 
experiment, it was not possible to resolve where the flame 
was relative to the burnout point, so the location of the 
flame in fig. 16b is a guess. The predicted spread rate of 
0.38 cm/s was below the experimental value of 0.55 cm/s*. 

Although the experimental data is not existent, 
various solid phase dimensional profiles from the computa- 
tion are presented in fig. 17. 

Finally, spread rate data at a free stream velocity of 
approximately 5 cm/s and at different oxygen concentrations 
are shown in fig. 18. The model is compared to the experi- 
ment. In the experiment, the velocity of both the flame 
base and the flame tip is reported, as shown in fig. 18, 
just before the end of the drop. Thus, the flames were 

*The experimental spread rate actually is the mean of 
the flame tip and flame base spread rates. They were close 
(0.53 and 0.57 cm/s, respectively) but the flame was still 
shrinking slightly as it approached steady state. The 
shrinking flame suggests the spread rate would slow down, 
given more time. 
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Figure 1 7. Dimensional profiles. 4o= 5 cm/s, Xo-1 5% 
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reported as either getting shorter or longer at the end of 
the drop. Reasonably good agreement is obtained, assuming 
the flames will eventually reach steady state. One reason 
for the difference in the quench limits may be due to the 
short drop time, as a flame may need more time to go out. 

An extensive comparison of all of the experimental 
data at this time is not practical since, as mentioned, 
many of the tests were still transient. The above data are 
presented to show that the model produces qualitatively the 
correct flame shape and spread rate trends, even though gas 
phase radiation, detailed flame chemistry, and a complicat- 
ed pyrolysis model were not considered. 

4.3. Parametric Comparison of Theoretical Results 

The two parameters varied in this work are free stream 
velocity and oxygen concentration. Several comparisons 
will be presented. 

4.3.1. Global Results: Spread Rates and Maximum Temperature 

In fig. 19, the computed spread rates and maximum gas 
phase temperatures are presented for all of the computed 
cases. The data are plotted against the reference veloci- 
ty, which is the relative gas velocity with respect to the 
solid burnout point (or flame) and is the quantity used in 
nondimensionalizing the equations. The spread rate 
increases approximately linearly with either flow velocity 
or oxygen percentage. Linear regression is used to get the 
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following relation, specific for the fuel used in this 
work: 

Vp = 1.12X 02 Ue.- 0.073IL + 0.194X O2 - 0.121 (84) 

In eq. (84), U*, an d V F are expressed in cm/s. In fig. 19, 
the dotted lines are flammability limits beyond which a 
solution does not exist. (The flammability limits will be 
discussed in a §4.4.) The spread rate trends and magni- 
tudes are qualitatively similar to those obtained in the 
experiment. 30 Another study 32 examined concurrent flow 
flame spread over a thin fuel in a normal gravity horizon- 
tal wind tunnel at higher speed (>35 cm/s) . After igni- 
tion, the flames would grow and eventually reach a steady 
state. Consistent with the theoretical results shown in 
fig. 19a, in ref. 32, the flame spread rate increased 
monotonically with forced velocity, in the velocity range 
less than 1 m/s. At high velocity (>1 m/s), however, the 
spread rate became constant as the flow velocity was 
increased. 

In fig. 19b, the maximum flame temperature, normalized 
by adiabatic flame temperature, is plotted. Near the 
quench limit (at low velocities) the maximum temperature 
drops off dramatically. This is a result of the increased 
relative importance of heat loss for these small flames. 
The rate of heat loss becomes a significant percentage of 
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the combustion heat release rate. The resulting drop in 
flame temperature corresponds to the decrease in spread 
rate, as the heat flux to the fuel is reduced. 

4.3.2. Parametric Comparison of Flame Structure 

In order to demonstrate the effect of oxygen and flow 
velocity on flame structure, two sets of cases will be 
presented, one at constant oxygen level and the other at 
constant velocity. 

In fig. 20, reactivity contours are shown in non- 
dimensional coordinates. The oxygen concentration is held 
at 15% and the flow velocity is varied from 2.7 to 15.5 
cm/s. The flame becomes longer as velocity increases. The 
flame strength also increases. This is evident in the 
larger size of the highest reactivity zone. Also, at 
higher velocity, the reactivity contours tend to bend back 
toward the fuel downstream, indicating increased flame 
strength . * 

The visible flame can be defined as the region 
enclosed by a reactivity contour (10 -4 g/cm 3 /s is assumed 
in this work) . For the simple reaction scheme in this 
model, the rate of fuel reaction is the best measure of how 
the flame appears. Based on these contours, then, the 


*The flame strength can be quantified as the area of 
the reactivity contours. This is a measure of the total 
heat release. 







83 


flame aspect ratio (i.e., the length to width ratio) 
changes dramatically with velocity. Additionally, the 
maximum reactivity is always located very near the burnout 
point, which is the most upstream part of the flame and 
sees the oxygen first . 

In fig. 21, the same flames are presented in dimen- 
sional coordinates. This plot demonstrates that flame 
length is greatly affected by velocity. In ref. 32, the 
flame length is reported to decrease with an increase in 
flow velocity (when is between 0.35 and 1 m/s) before 
leveling off at high velocity (greater than 1 m/s) . While 
contrary to the predicted results, this may be an effect 
evident only at higher velocities. 

The effect of different oxygen concentrations at a 
fixed free stream velocity is shown in fig. 22. The flame 
length decreases slightly with oxygen. Again, the flames 
at higher oxygen are stronger based on the fact that they 
are bigger and tend to curve back toward the fuel down- 
stream. In fig. 23, the flames are shown in dimensional 
coordinates. The difference between fig. 22 and fig. 23 is 
small because the relative velocity does not change to much 
for the different cases. Hence, the thermal lengths are 
approximately equal . 

4.3.3. Flame Stand-off Distance and Thickness 

When the flame thickness and stand-off distance are 
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Figure 21 . Dimensional fuel reactivity contours at 1 5% O 2 . Weakest 
(outermost) contour is 1 0~ 6 g/ cm 3 / s. 

A factor of 1 0 separates contours. 
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defined as shown in fig. 24, a comparison can be made. 
Fig. 24 presents both dimensional and nondimensional 
results. Looking at part (a) first, the dimensional flame 
stand-off distance and thickness increase with decreasing 
flow velocity except very near the quenching extinction 
limit. This is explained as follows. Initially, as the 
flow velocity is reduced, the flame moves away from the 
fuel and thickens, as it tries to move closer to the region 
of fresh oxidizer. As the flame moves away, the heat flux 
to the solid diminishes, leading to a decrease in the rate 
of pyrolysis. This makes the flame shorter and increases 
the relative importance of heat loss, hence cooling the 
flame. When the flame dimensions become the same order as 
the thermal length, two-dimensional heat conduction loss 
becomes important, further decreasing the flame tempera- 
ture. When this heat conduction loss becomes important, 
the flame forsakes the search for more oxygen and actually 
moves back toward the fuel in an attempt to maximize the 
ratio of heat flux to the solid to heat loss to the 
environment . 

In fig. 24b, the dimensional flame stand-off distance 
increases monotonically with oxygen percentage. This is 
due to the increase in flame temperature with oxygen 

*In a three-dimensional problem, the conductive loss 
in all three dimensions would become important. 
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percentage. At higher oxygen, the flame simply does not 
have to get as close to the fuel because it is hotter. The 
burning rate is also higher, pushing the flame farther from 
the fuel. The flame is thinner at higher oxygen percentage 
because the mass diffusion driving potential is higher. 
Ultimately, the flame goes out as the oxygen percentage is 
reduced through the finite-rate chemical reaction term. 
The extinction mechanisms will be described in §4.4. 

4.3.4. Solid Phase Parametric Results 

The characteristics of the solid fuel are compared in 
fig. 25 which shows the effect of varying forced flow 
velocity on solid temperature and incident heat flux. In 
part (a) , the surface temperature plots indicate that the 
maximum always occurs at the burnout point, as expected. 
Furthermore, the maximum increases as the flow velocity 
increases. This is due to the fact that the flame is 
pressed closer to the fuel at high velocities, increasing 
the heat flux and thus increasing the temperature. The 
heat flux profiles are shown in part (b) . The increase in 
maximum heat flux due to increasing flow velocity is shown. 
Note the sharp bend in the curves . The point where the 
bend occurs indicates the end of the pyrolysis region. In 
general, as flow velocity is decreased, the heat flux, 
solid temperature, and pyrolysis length decrease. For the 
lowest-energy flames in this work (i.e. low fuel consump- 




Figure 25. Solid temperature and heat flux profiles at 1 5% O 2 
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tion rate or spread rate) , the pyrolysis region begins to 
approach point source behavior, as quenching extinction is 

neared (described in §4.4) . 

The heat flux in all cases varies approximately as x , 
where n— 1/2 in the pyrolysis region, and n— 1 in the 
preheat region.* Considering a non-reacting boundary layer 
with heat transfer, for an isothermal plate, q„~ * 1/2 • This 
is approximately given in the pyrolysis region. In 
general, if the appropriate temperature difference in the 
boundary layer is specified as (T g -T w ) ~x m , then q w ~x 
Thus, if m=-l/2, the heat flux to the solid in the preheat 

region for this problem is reasonable. 

In fig. 26, similar results are plotted to demonstrate 
the effect of varying oxygen percentage at fixed free 
stream velocity. Again, the pyrolysis length is clearly 
shown by the bend in the heat flux curves. The solid 
temperature, heat flux, and pyrolysis length, decrease with 
decreasing oxygen percentage. (In fig. 26b, the curve at 
21% 0 2 has a small kink right at the elliptic/parabolic 
interface, x=40 cm. This is a manifestation of the 
difference between the equations in the two domains.) 
4 . 3 . 5 . Definition of Various Length Scales 

The definition of pyrolysis length is somewhat 

*The actual exponents are closer to n=-0.43 in the 
pyrolysis region and n=-l.l in the preheat region. 
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arbitrary. In an experiment, the pyrolysis front may be 
determined as the point where blackening of the paper first 
occurs. 32 In a theoretical model, the pyrolysis front is 
often described as the point where a certain percentage of 
fuel remains. In this work, the pyrolysis front was chosen 
as the point where 99.9% of the fuel remains. It turns out 
that for all cases, this specification corresponds very 
closely to the bend in the heat flux curves (seen in fig. 
26, for example) . 

Another length scale of interest in the solid phase is 
the preheat length, which is the distance required for the 
fuel to heat up to the temperature at which pyrolysis first 
occurs. in this work, the definition of preheat length is 
based on the heat flux to the solid, as follows. The heat 
flux at the end of pyrolysis is found. The preheat length 
is defined as the point from the end of pyrolysis to the 
point where the heat flux has dropped by a factor of ten. 

Finally, the flame length is defined as the length of 
the visible flame, which was defined earlier as the region 
bounded by the fuel reactivity contour of 10" 4 g/cm 3 /s. 

The results are summarized in fig. 27. The pyrolysis, 
flame and preheat lengths all vary linearly with flow 
velocity and oxygen concentration. 

4.4. Extinction Boundary 

One of the main goals of this work was to determine 




Figure 27. Flame ( l F ), pyr 
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the extinction characteristics of this flame spread 
problem. Work done in the past has demonstrated that 
flames in very low speed flows will eventually go out since 
they become weakened due to heat loss (see ref. 20, for 
example) . This mode of extinction is called quenching. 
Another more familiar mode of extinction, called blowoff, 
occurs when the flow velocity becomes too high. The 
chemical reactions are too slow compared to the rate of 
incoming oxidizer, so the flame cannot be stabilized and is 
blown off. 

At normal gravity, there is a limit to how small the 
velocity of oxidizer into the flame zone can be. This is 
due to buoyancy, which sets up flow velocities on the order 
of several tens of centimeters per second. At first, our 
experience with blowoff in normal gravity suggests that if 
we first removed buoyancy, then reduced the flow velocities 
to arbitrarily small values, the flame would always exist. 
This is because the chemical reactions have relatively more 
and more time to proceed at the lower flow speeds. What 
actually happens is that the reduced flow of oxidizer 
diminishes the overall power of the flame in that heat 
losses (due to radiative loss, for example) drop the flame 
temperature below the point that the reaction can be 
sustained. This is the quenching extinction mechanism. 

In fig. 28, the extinction boundary is presented. The 
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40 circles represent flame solutions, and the 9 X's show 
the extinction points. The two extinction branches, 
quenching and blowoff, are shown. In obtaining the 
extinction points, care was taken to make small changes in 
velocity or mole fraction. In a steady model such as this, 
there is always a possibility that the extinction point is 
actually a non-convergence point. By making small changes 
in the parameters, however, it is possible to get arbi- 
trarily close to the actual extinction boundary. 

In fig. 28, a logarithmic scale is used to plot free 
stream velocity so that detail at low velocity is evident. 
A plot using a regular scale is given in fig. 29. The 
quenching and blowoff branches have dramatically different 
slopes on this scale. 

4.4.1. Differences Between Blowoff and Quenching 

The blowoff and quenching branches have different 
characteristics. While the equations used are steady, some 
indication of an unsteady time response is obtained by 
looking at the results during iteration.* The blowoff 
phenomenon is indicated in fig. 30. The qualitative, 
temporal response is correct . Here, the top view shows 
reactivity contours for a converged flame near the blowoff 
boundary. As the blowoff boundary is crossed (by increas- 

*Underrelaxation is employed, so the iteration number 
is not exactly proportional to time. 




Figure 29. Extinction boundary plotted in regular coordinates. 
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ng velocity or decreasing oxygen concentration) the 
sequence of plots show that the flame is unable to stabi- 
lize and is blown downstream. In the second view, the 
flame starts to lift off at x-0. This flame is actually 
stronger (based on increased area of reactivity contours) 
than the one above it, due to the penetration of oxygen 
beneath the flame. The remaining views show the flame 

getting weaker and moving farther downstream. Eventually, 
the flame is blown out. 

In contrast, quenching extinction is observed to be 

very different. If the quenching boundary is crossed (by 

decreasing velocity or oxygen concentration) the flame 

always remains firmly anchored to the burnout point (x=0) 

but merely gets weaker and smaller until it disappears 
entirely . 


There are other clues Indicating the different nature 
of the two extinction modes. Comparing a near quenching 
and near blowoff flame, fig. 31 depicts the isotherm 
structure and maximum temperature. At 2.7 cm/s, the flame 


has its maximum temperature 
burnout point. At 15.5 cm/s, 
blown downstream somewhat, 
dramatic in nondimensional 


slightly upstream of the 
the maximum temperature is 
The effect is even more 


coordinates: the maximum occurs 



f m n i 


(a) X 0 =1 57,, I4p=2.7 cm/s 



(b) X 0 =1 57, LU=1 5.5 cm/s 

Figure 31 . Isotherms for a near quench point (a) 
and near blowoff point (b) at 1 5% O 2 . 
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over five thermal lengths downstream of the burnout point.' 

To get a better understanding of the overall flamma- 
bility boundary, fig. 32 presents maximum temperature and 
spread rate for points on the boundary. The inset shows 
the boundary and the points used. The maximum temperatures 
on the quenching branch are quite low. Furthermore, it 
appears that a near-limit flame will tolerate a larger 
reduction in temperature at higher oxygen concentrations. 
At blowoff, the maximum temperature is a healthy percentage 
of the adiabatic flame temperature. Ouenching extinction, 
on the other hand, is ultimately due to weakened chemical 
kinetics through flame temperature reduction.' 

In both blowoff and quenching extinction, the flame 
eventually goes out due to weakened chemical kinetics. In 
the nondimens ional expression for reaction rate, given by 
eq. (33), the Damkohler number. Da, appears. As described 
earlier, Da is the ratio of the characteristic flow time to 
the characteristic chemical reaction time, in one thermal 
length. At blowoff, the reaction does not have time to 
proceed because of the high flow velocity. in this case, 

very near the' burnout point’. ^ maxilnunl reactivity occurs 

quenching 0 extinction, k p?ldictin P g la the an bo im d° rtant r0le in 
conf ident^that we^e' 

adequate to capture the essential “tu“7 ofThe p'roWem? 


ction Boundary 
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Da becomes too small. At quenching. Da is large, but the 
reduced flame temperature due to heat loss slows the 
reaction rate. Thus, while the mechanisms are quite 

^ eren *'' the flame goes out in either case as a result of 
the finite— rate reaction. 

There are other characteristics of the near blowoff 
(high velocity) flames. In fig. 33, two effects are shown. 
First, look at the reactivity contours near the fuel in the 
region from x=0 to x=3 . These contours are distorted due 
to the large flux of pyrolysis products. The concentration 
is so large that the reaction rate is noticeably enhanced. 
Second, follow the Y F =0.001 contour from small to large x- 
values. At first, the line indicates that some fuel has 
escaped the reaction zone since it lies outside the 
indicated reactivity contours. Eventually, the fuel is 
blown back into the flame where it is consumed. This 

effect is also seen to a lesser extent in the Y F =0.01 
contour . 

4.4.2. Incomplete Combustion 

Another interesting aspect relates to the flux of fuel 
in the x-direction as a function of x. In fig. 34, fuel 
flux (normalized by the total burning rate) is plotted for 
all cases at 15% 0 2 . Looking at an individual profile, the 
flux of fuel is incremented by the flow of pyrolysis gases 
from the solid. At the same time, the flame consumes the 
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fuel and eventually leads to a drop in fuel flux. Finally, 
the reaction ceases and an amount of fuel escapes from the 
flame unreacted. This "escaped fuel" or (incomplete 
combustion) is plotted in the inset. As quenching is 
approached, the amount of escaped fuel increases dramati- 
cally. There are three reasons the escaped fuel increases 
near quenching. First, the relatively cool fuel plate 

quenches the gas phase reaction near the plate. Second, 
the radiative heat loss from the solid fuel leads to 
reduced gas phase temperature. Third, the short flames 
simply permit an increased amount of fuel to escape. This 
is based on simple geometric considerations, namely, the 
flame length shrinks to about the same order as the flame 
stand-off distance. 

As a final word, blowoff at higher oxygen levels was 
not searched for in this work. There are three reasons. 
First, this work was mainly interested in low speed 
quenching effects. Second, flames at high velocity eventu- 
ally become turbulent, making this laminar model inade- 
quate. Third, in some cases, a near-blowoff flame can 
become side stabilized. The flame retreats downstream of 
the burnout point, where the velocity is lower. After 
burning through the fuel at this location, it retreats 
again, and repeats the cycle. Clearly, this steady model 
is insufficient in this case. 



Conclusions 


. A numerical model of concurrent flow flame spread over 
a thin fuel has been formulated. The effect of oxygen 
concentration and forced flow velocity on flame structure 
and extinction limits in a zero gravity environment has 
been found. 

The governing equations included the steady, laminar, 
two-dimensional, fully elliptic conservation equations for 
momentum, energy, and species. To capture the long flame 
and preheat zone without using an elliptic formulation 
everywhere, parabolic equations were solved at a point 
sufficiently downstream from the leading edge region. The 
parabolic solution provided the boundary condition for the 
elliptic equations, minimizing the effect of the transition 
between the two sets of equations. 

To get a steady solution, the coordinate system was 
chosen to move with the flame. Thus, the flame was 
stationary and the fuel fed in at the rate of flame spread. 

The flame structure has been presented in detail. 
Contours of temperature, fuel reactivity rate, and equiva- 
lence ratio have been compared. The solution of the 
Navier-Stokes equations allowed the velocity and pressure 
fields to be visualized. Streamlines showed the deflection 
of the flow due to the presence of the flame and fuel, as 
well as the flow of fuel vapor emanating from the solid due 
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to pyrolysis. Solution of the solid phase equations gave 
the pyrolysis length, the incident heat flux distribution, 
and of course the solid temperature. 

A comparison was made with experiment, which reported 
visible flame shape and spread rate. Since the experimen- 
tal data was in many cases still approaching steady state, 
only one nearly steady flame was compared. Using a contour 
of constant fuel reactivity to represent the visible flame, 
the flame shape was predicted quite accurately. The spread 
rate was under-predicted by about 20% for the case exam 
ined- No attempt was made to adjust the parameter values 
in the model to more closely match the experiment, but 

results agreed qualitatively. 

An extinction boundary, with flow velocity and oxygen 
concentration as coordinates, was presented. It consisted 
of two branches, a blowoff branch and a quenching branch. 
This type of extinction boundary has been observed experi- 
mentally in opposed flow flame spread tests. 

The flame went out in different ways depending on 
which branch of the extinction boundary was crossed. In 
blowoff, the chemical reaction did not have time to proceed 
because of the high rate of flow, and the flame was subse- 
quently blown downstream. The flame temperature was high 
enough to allow the chemical reaction to proceed, however, 
it was found that the maximum temperature zone in the flame 
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was pushed downstream of the burnout point and higher off 
the fuel surface as the blowoff branch was approached. 
Therefore, the delayed chemical reaction effectively 
reduced the heat flux back to the fuel in the burnout 
region. Thus, the required amount of fuel vapors needed to 
stabilize the flame could not be generated, and the flame 
was eventually blown downstream. 

As the quenching branch was approached, on the other 
hand, the flame was always firmly stabilized at the burnout 
point . Furthermore, the maximum temperature in the flame 
always occurred very near the burnout point. Thus, the 
chemistry had plenty of time to proceed. However, the 
maximum flame temperature decreased steadily as the 
quenching branch was neared. This was due to the increased 
importance of heat loss (radiation from the solid fuel) as 
the visible flame size and power output (equivalently, 
rate) were reduced. The flame was quenched when the 
temperature became too low to sustain combustion. Addi- 
tionally, flames near the quench limit permitted an 
increased percentage of fuel vapors to escape unreacted, an 
indication of the relatively low flame temperatures and 
short flames. 

By looking at cases which go to extinction during 
iteration, an indication of the unsteady time response of 
the flame could be obtained, even though the model was for 
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steady solution. These results showed that blowoff was 
characterized by the flame being de-stabilized and pushed 
downstream, while quenching was due to the gas phase flame 

becoming cooler and cooler. 

The theoretical results are compared. At constant 
oxygen concentration, the affect of flow velocity from 
quenching to blowoff was presented. At constant free 
stream velocity, the effect of different oxygen levels, 
from 21% (air) all the way down to the extinction limit, 
was shown. Flame length, pyrolysis length, preheat length, 
and spread rate increased approximately linearly with both 
free stream velocity and oxygen concentration. The maximum 
temperature in the gas phase decreased as flow velocity was 
reduced, and dropped off steeply as the quenching branch 


was approached. 

From one case to another, the solid temperature and 
incident heat flux profiles had similar shapes. Pyrolysis 
and preheat regions were identifiable. The biggest differ- 
ence was that the maximum temperature and incident heat 
flux (which always occurred right at the burnout point) in 


creased with flow velocity and oxygen percentage. 

This model provides a powerful tool in examining flame 
spread and extinction characteristics. It is a solid base 
upon which further modifications can be made. Eventually, 
it will be generalized to consider unsteady flame spread 
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over more common (thicker) fuels. 

Other issues worth considering include the following. 
Gas phase radiation should be included to get a better 
gauge of its effect on flame spread, structure, and 
extinction (as a heat loss mechanism). An improved gas 
phase chemical reaction formulation is needed to model the 
low temperatures near quenching, where intermediate species 
could be important. The presence of soot will have an 
influence on both the gas phase chemistry and radiation. 
Finally, the solid fuel model should be improved to more 
realistically predict the rate and type of pyrolysis 
products being generated. 
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Appendix 


The FORTRAN computer program used to generate these 
results follows. In addition, the important parameter 


values (read in unit 10 in the program) are listed below. 


Ill 
4 4 4 

10 3 

69 49 26 

150 1 1 

7 . OOOOOE-Ol 
6 . 50000E-01 
1 . OOOOOE+OO 
1.18500E+00 
2 . 75000E-04 
3 . OOOOOE+02 
9 . 09000E-01 

70 20 65 

0 1 65 

30 5 1 


11111 
4 4 0 1 

4 0 

1.1000 1.1000 0.2000 
- 0.0001 - 0.0002 - 0.0001 
9. 00000E-01 6 . 50000E-01 

6 . 50000E-01 6 . 50000E-01 

4 . 16667E+00 9.25485E+08 

2 . 10000E-01 2. OOOOOE+OO 

2 . 63000E-01 3. OOOOOE-Ol 

3. 00000E+02 -1 . 80000E+02 
1. OOOOOE+OO 3 . OOOOOE+Ol 
1.1000 3.0000 30.0000 

0.0000 2.0000 


0.2000 3.0000 3.0000 

- 0.0001 - 1.0000 - 0.0012 


6 . 50000E-01 
3 . 80000E+07 
4 . 04000E+01 
2 . 38254E-01 
3 . 30000E-01 
3 . 80000E-03 
1. 00000E+00 


5 . 03000E+01 
4 . 53000E+01 
1. 00000E+00 
1 . 35600E-12 
2 . 13000E+00 
0. 00000E+00 
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PROGRAM FLAME 

COMMON/L1/DX (0:93) , DY (0 : 50) , NX, NY,NLE 
COMMON/L2/U <0:93, 0:50) ,V (0:93, 0:50) ,T (0:93, 0:50, 1:3) , 
R(0:93,0:50) ,P (0:93, 0:50) ,PPR (0:93, 0:50) 
COMMON/L3/UHAT(0: 93, 0:50) , DU (0 : 93, 0 : 50) ,USTAR(0: 93, 0:50), 

UAP (93,50) ,UAE (93,50) ,UAW(93, 50) , 

UAN (93, 50) ,UAS (93, 50) ,UB (93, 50) , GRASH, RINFNON 
COMMON/L4/VHAT (0 : 93, 0 : 50) ,DV (0 : 93, 0 : 50) , VSTAR(0 : 93, 0 : 50) , 

VAP (93, 50) , VAE (93, 50) ,VAW(93,50) , 

VAN (93,50) , VAS (93, 50) , VB (93, 50) 

COMMON/L5/AP (93, 50) ,AE(93,50) , AW (93, 50) , AN (93, 50) , 

AS (93, 50) , SC (93,50) ,SP (93,50), B (93, 50) 

COMMON/ L6/AA (0 : 93) , BB (0 : 93) ,CC (0 : 93) ,DD (0 : 93) , XX (0 : 93) 

COMMON /L7 /RSTAR, RS, CS, CP, SIG, TO, TL, XL, ES, TAU, ASTAR, TSTAR 
COMMON/L8/DXP (0:500) ,DYP (0:100) ,UPAR(2, 0:100) , VPAR (2, 0:100), 

TPAR (2, 0:100,3) ,TS(0:570,2) ,USPAR(0 : 570) ,VSPAR(0 : 570) 
COMMON/ L9/ IBC, UBC (0:50) , VBC (0 : 50) , TBC (0 : 50, 3) 

REAL RX, RY, DXMIN, DYMIN, DXMAX, DYMAX, TAMB, DA, Q, E, PR, FO, YOAMB, 

— ERM, ERP , ERU, ERV, ERPP, ERF, PDAM, UDAM,VDAM, FDAM ( 3 ) ,VEL,TOVI (0:570) , 
-POLD (0 : 93, 0 : 50) , PH (2) ,MAINPDIF, THICK (0 : 570) , DAD, VE0, VF, EPS, 
-ASW, VFINI , QO, QN, QOL, QNL, GLEVEL, GRASH, RINFNON, UBUOY 

C GET GENERAL DATA 

READ (10, 910) IPC, IUC, IVC, IPPC, ITC, IYFC, IYOC, IOUV 

READ (10,910) IBP, IBU, IBV, IBPP, IBF, IUSTVST, IDA 

READ (10, 910) ITSSMOO, ISBC, JSK, KSOL, ISKPP 

READ (10, 900) NX, NY, NLE,RX,RY, DXMIN, DYMIN, DXMAX, DYMAX 

READ (10, 900) NALL, IREADOLD, NSCH, ERM, ERP, ERU, ERV, ERPP, ERF 

READ (10, 905) PR, PDAM, UDAM, VDAM 

READ (10, 905) FDAM ( 1 ) , FDAM ( 2 ) , FDAM ( 3 ) ,ASW,ES 

READ (10, 905) TAMB, TSTAR, DAD, Q,E 

READ (10, 905) FO, XOAMB, VE0 , VF, EPS 

READ (10, 905) RSTAR, RS,CS, CP, SIG 

READ (10, 905) TO, TL, XL, TAU, ASTAR 

READ (10, 905) QO,QN, QOL, QNL, GLEVEL 

READ (10, 900) IBEG, JTOP, I RIGHT, RDXP, DXPMIN, DXPMAX 

READ (10, 900) ISKIP, ISR, IPARN, VFL1 , VFU1 

READ (10, 900) MMAX, KPAR, IBC 

YOAMB=l. / (1 .+ (1 . /XOAMB— 1 . ) *0.8754513) 


C GENERATE GRID 

CALL GRID (RX,RY, DXMIN, DYMIN, DXMAX, DYMAX) 
C GENERATE INITIAL FIELD 


CALL READIN (TAMB, YOAMB, TSTAR, IREADOLD, IUSTVST, ISR, IBEG, IRIGHT) 

c TAILOR INITIAL SURFACE TEMPERATURE DATA 

IF (ISR.EQ.l) THEN 

IF (IRIGHT.GT. IPARN) IRIGHT=IPARN 
IF (IRIGHT.LT. IPARN) THEN 

DO 3 I=IBEG+IRIGHT+1 , IBEG+IPARN 
TS (I, 1) =TS (IBEG+IRIGHT, 1) 

3 TS (I, 2) =TS (IBEG+IRIGHT, 2) 

IRIGHT=IPARN 

ENDIF 

DO 4 J=NLE+1, IBEG+IRIGHT 

C :*************************************************************.* 

C** NOTE!!! THE FOLLOWING LINE ASSUMES EQU. OF STATE ** 

C* *********************** RHO=TSTAR/T *********************** 

4 VSPAR ( J) =EXP (-ES/TS ( J, 1 ) ) *ASW*RS/TSTAR*TS ( J, 1 ) /RSTAR/ 

(VEO+UBUOY-VF) 


ENDIF 
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motf.- it is important that ibeg+iright>- nx+1- - 

1 main routine 

NF=0 

NG=0 

IPAR=0 

IOKO=0 

VFRI=VF 

QORI=QO 

» »»»« ° F STOTE 

RI NFNON=T S TAR / T AMB 
VEL=VE 0 +UBUOY-VF 
DA=DAD/VEL/VEL 
GRASH= (UBUOY /VEL) **3. 

CHERPP=0 . 04 

IF (ERPP . EQ. -1 . ) CHERPP=0 . 

CHQO=0 . 

IF (QN.EQ.l.) CHQO=l . 

DO 200 1=1, NALL 

VEL=VEO+UBUOY-VF 
DA=DAD/VEL/VEL 
GRASH= (UBUOY/VEL) **3. 

IF (IOUV.EQ.0) GO TO 5 

CALL OLDS (POLD, P , NX, NY) 

CALCULATE UHATS AND VHATS 

CALL UHATS (PR, TSTAR) 

CALL VHATS (PR, TSTAR) 

SOLVE FOR THE PRESSURE FIELD- - 

5 DO 10 K=1 , IPC 

10 CALL PRESSURE (ERP,PDAM, IBP) 

CALCULATE USTAR AND VSTAR 

DO 20 K=1 , IUC 

20 CALL USTARS (ERU,UDAM, IBU, PR, TSTAR) 

DO 30 K=1 , IVC 

30 CALL VSTARS (ERV,VDAM, IBV, PR, TSTAR) 

SOLVE THE P' EQUATION 

DO 40 K=1 , IPPC 

CALL PPRIME (ERPP, IBPP, ISKPP) 

CORRECT VELOCITY ~ 


40 CALL CORVEL (PR, TSTAR) „„ nTrI „ 

COMPUTE OXYGEN CONCENTRATION FIELD- 


DO 50 K=1 , IYOC _ 

50 CALL PHI (TAMB, DA, Q,E,FO, ERF, FDAM, 3, TSTAR, IBF) 

c COMPUTE FUEL CONCENTRATION FIELD 

DO 60 K=1 , IYFC 

60 CALL PHI (TAMB, DA, Q, E, FO, ERF, FDAM, 2 , TSTAR, IBF) 

COMPUTE TEMPERATURE, DENSITY - 

DO 70 K=1 , ITC 

CALL PHI (TAMB, DA, Q,E,FO, ERF, FDAM, 1, TSTAR, IBF) 

u -_ I------- DEN S < TS TA R) , SET up GRID and INITIAL PARABOLIC VALUES 

IF ( (IOKO.NE. I/KPAR) .OR.I.EQ.l) THEN 
IF (IBEG.GT.NX+1.OR.IRIGHT.LE.0) THEN 
DO 80 I4=NLE+1 , NX+1 
TS (14, 1) =T (14,0,1) 

80 TS (14 , 2 ) =T (14,1, 1) 

E Sll PARGRID (NF, IBEG, JTOP, IRIGHT, RDXP,DXPMIN,DXPMAX) 
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83 


87 

90 


100 


102 


C) ) ) ) 


IF (NG.EQ.l) GO TO 83 
NG=1 

Tf L ! , T S rI S mT <I ? EG ' JTOP ' IRIGHT, ES,ASW,RS,TSTAR,RSTAR,VF,VEL, ISR) 

±r ( I. EQ. NALL) IPAR=1 
ENDIF 
ENDIF 

IF ( (IOKO.NE. I/KPAR) .OR. I .EQ. 1) THEN 

COMPUTE SOLID TEMPERATURE AND THICKNESS 

IF (NSCH.EQ. 1) THEN 
DO 87 IK=NLE+1, I BEG 
TS (IK, 2) =T {IK, 1,1) 

DO 90 IK=NLE+1, IBEG+IRIGHT 
TOVI(IK)=TS(IK,l) 

VFINI=VF 

ZZZ=0. 

VFL=VFL1 

VFU=VFU1 

CALL SOLID (EPS, VFL,VFU, VF, VE0,UBUOY, 

nn inn t mt ^i IC t4 ZZZ ' TX ' Q0 ' QN ' ASW ' IBEG ' IRIGHT, MMAX, KSOL) 

DO 100 J=NLE+1, IBEG+IRIGHT 

TS (J, 1) = (QNL*TS (J, 1) +QOL*TOVI (J) ) / (QNL+QOL) 

IF (ITSSMOO.NE. 0) THEN 

DO 102 J=NLE+2, IBEG+IRIGHT 

IF (TS ( J, 1 ) . GT . TS ( J— 1, 1 ) ) TS (J, 1)=TS (J— 1, 1) 

ENDIF 

) ) ) ) COMPUTE SPREAD RATE BY INTEGRATING MASS FLUX ( ( ( ( 
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XINTVF=0 . 

DO 103 L=NLE+1 , IBEG+IRIGHT— 1 
IF (L.LT.IBEG) DDX=DX(L) 

IF (L.EQ. (NLE+1) ) DDX=DX(L)/2. 

IF (L.GE.IBEG) DDX=DXP (L-IBEG+1) 

C_ 7«o N 0TE: ™ E integra L BELOW ASSUMES SQRT (1+ (DH/DX) **2) = 1 

103 XINTVF=XINTVF+DDX*EXP (-ES/TS (L, 1) ) 

XINTVF=ASW / TAU* ASTAR/ (VEO+UBUOY-VF) *XINTVF 

QO=QO*XINTVF/VF*CHQO+QO* (l.-CHQO) 

WRITE (56, 104) VF, XINTVF 
VF= (QNL*VF+QOL*VFINI) / (QNL+QOL) 

FORMAT (IP, 2E12 . 5) 

C CALCULATE BLOWING VELOCITY 

DO 105 J=NLE+1, IBEG+IRIGHT 

C************************************************************** 

^!***!*** NOTE!!! THE FOLLOWING LINE ASSUMES EQU. OF STATE ** 

° *********************** RHO=TSTAR/T *********************** 
105 VSPAR ( J) -EXP (-ES/TS (J, 1) ) *ASW*RS/TSTAR*TS (J, 1) /RSTAR/ 

- (VEO+UBUOY-VF) ' 

VSPAR (NLE+1 ) =VSPAR (NLE+1 ) /2 . 

C SET B.C. FOR SOLID FUEL 

IF (ISBC.EQ. 0) THEN 

TS (IBEG+IRIGHT, 1) =TS (IBEG+IRIGHT— 1, 1) 

ELSE 

TS (IBEG+IRIGHT, 1 ) =1 . 

ENDIF 

IF (ZZZ.EQ.l.) GOTO 210 
ENDIF 

IF (JSK.EQ.0) GO TO 185 

C SOLVE PARABOLIC REGION 

IF (IBEG. LE.NX+1 .AND. IRIGHT. GT. 0) 

CALL PARCA (IBEG, JTOP, IRIGHT, PR, TSTAR,GR, RINFNON, DA, Q, E, FO, 
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C 


IPAR, IDA, VEL, RSTAR, ASTAR) 

CORRECT BOUNDARY CONDITIONS FOR U-VELOCITY 

VMEAN= (VSPARai) *DX{li) +VSPAR (II+l ) *DX (II+l) ) / (DX(II) +DX (II + l) ) 
DHDX= (THICK (II+l) -THICK (II) ) / (DX (II+l) +DX(II) ) *2. 

DHQ=DHDX/SQRT (1 . +DHDX*DHDX) 

110 USPAR (II ) =— VF/ (VEO+UBUOY-VF) — DHQ*VMEAN 
DO 120 II=IBEG/ IBEG+IRIGHT— 1 
VME AN= V S P AR (II) 

IF ( (II.EQ.IBEG) -AND. (IBEG.EQ.NX+1) ) THEN 
DXL=DX (NX) /2 . 

DXU=DXP ( 1 ) 

GOTO 112 
ENDIF 

IF (II.EQ.IBEG) THEN 
DXL=DX (IBEG) / 2. 

DXU=DXP ( 1 ) 

ELSE 

DXL=DXP (II-IBEG) 

DXU=DXP (II— IBEG+1) 


ENDIF 

112 RREP=DXL/ (DXL+DXU) 

DHDX= (THICK (II+l) -THICK (II) ) /DXU*RREP+ 

- (THICK (II) -THICK (II-l) ) /DXL* (1 . — RREP) 

DHQ=DHDX/SQRT (1 .+DHDX*DHDX) 

120 USPAR ( II )=-VF/ (VEO+UBUOY-VF) -DHQ*VMEAN 
USPAR (IBEG+IRIGHT) =USPAR (IBEG+IRIGHT— 1) 


C 


C 


II=NLE 

VMEAN=VSPAR (II + l ) 

DHDX= (THICK (II+2) -THICK (II+l) ) / (DX ( II+2 ) +DX (II+l ) ) *2. 
DHQ=D HDX / SQRT (1 . +DHDX*DHDX) 

USPAR (II) =-VF/ (VEO+UBUOY-VF) — DHQ*VMEAN 

UPDATE ELLIPTIC SURFACE VARIABLES 

(NSCH.NE.l) GO TO 155 


IF 

DO 130 II=NLE+1 , NX+1 
T ( 1 1 , 0 , 1 ) =T S ( 1 1 , 1 ) 

130 V(II, 0)=VSPAR(II) 

DO 140 II=NLE+1, IBEG-1 
140 U(II, 0)=USPAR(II) 

U (NLE, 0 ) =USPAR (NLE) 


DO 150 II=IBEG, NX 

RR1=DX (II + l) / (DX (II) +DX (II + l) ) 

150 U(II, 0) =USPAR(II) *RR1+USPAR(II+1) * (1.-RR1) 
155 CONTINUE 
ENDIF 

I OKO= I / KP AR 

UPDATE LEFT AND TOP VELOCITIES — 


DO 170 11=0, NX 

170 U (II , NY+1 ) = (VE0-VF) / (VEO+UBUOY-VF) 
DO 180 JJ=0 , NY+1 

180 U(0, JJ)=(VE0-VF) / (VEO+UBUOY-VF) 


185 IF (ISKIP.NE.0) GO TO 200 

CALL FINDDIF (P , POLD, NX, NY , MAINPDIF P , MAIN IT . , IDP , JDP ) 

PH (1) =PH (2) 

PH (2 ) =ABS (MAINPDIF) 

IF (PH(2) .GT.PH(l) ) THEN 
ICOUNT=ICOUNT+l 
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ELSE 


ICOUNT=0 

ENDIF 

IF ( (ICOUNT^GE^IO) .AND. (PH(1) .GT.0.1) ) GO TO 210 

FIND LARGEST, SMALLEST P PRIME 

PPMAX*— 100000. 

PPMIN=+100000 . 

DO 190 12=1, NX 
DO 190 J2=1,NY 

IF (PPR(I2, J2) .GT.PPMAX) THEN 
PPMAX=PPR (12 , J2 ) 

IDMA=I2 


JDMA=J2 

ENDIF 

IF (PPR(I2, J2) .LT.PPMIN) THEN 
PPMIN=PPR (12 , J2 ) 

IDMI=I2 
JDMI=J2 
ENDIF 
190 CONTINUE 

WRITE (22,915) PPMAX, IDMA, JDMA, PPMIN, IDMI , JDMI , 
IF ( (PPMAX-PPMIN) .LT.CHERPP) GO TO 210 
200 IF (ABS (MAINPDIF) .LT.ERM) GO TO 210 

WRITE DATA 

210 CONTINUE 


(PPMAX-PPMIN) 


CALL WRITE3 (T, 0, NX+1, 0, NY+1, 2) 

CALL WRITE2 (U, 0,NX, 0,NY+1, 2) 

CALL WRITE2 (V, 0,NX+1, 0,NY, 2) 

CALL WRITE2 (P, 1 , NX, 1 , NY, 2 ) 

CALL WRITE2 (PPR, 1, NX, 1 , NY, 2) 

CALL WRITE2 (USTAR, 0,NX, 0, NY+1, 2) 

CALL WRITE2 (VSTAR, 0, NX+1, 0, NY, 2) 

CALL WRITE 4 (TS, NLE+1 , IBEG+IRIGHT, 1,2,2) 

WRITE (13, 905) (UBC(J) , J=1,NY) 

WRITE(13, 905) (VBC(J) , J=1,NY) 

WRITE (13, 905) (TBC (J, 1) , J=1,NY) 

WRITE (13, 905) (TBC (J, 2) , J=l, NY) 

WRITE(13, 905) (TBC ( J, 3) , J=1 , NY) 

CALL WRITE 1 (THICK, TX, NLE+1 , IBEG+IRIGHT, 2) 

IF (IDA.EQ.l) THEN 

TO SAVE SPACE, DA DATA WRITTEN TO PPR 

PREF=DA*VEL*VEL*RSTAR/ASTAR 
DO 220 1=0, NX+1 
DO 220 J=0 , NY+1 

PPR (I, J) =PREF*R (I,J)*R(I,J)*T(I,J,2)*T(I,J,3) *EXP (-E/T(I,J,1) ) 
IF (PPR(I, J) .LT.l.E-78) PPR ( I, J) =0 . ' 

220 CONTINUE 

DO 230 1=0, NX+1 

230 WRITE(26, 905) (PPR (I, J) , J=0, NY+1) 

ENDIF 


CLOSE (UNIT=1 0 ) 

OPEN (UNIT=10) 

IF (ZZZ.EQ.l.) VF=VFRI 
IF (ZZZ.EQ.l.) QO=QORI 

REWRITE GENERAL DATA 

WRITE (10, 910) IPC, IUC, IVC, IPPC, ITC, IYFC, IYOC, IOUV 
WRITE (10, 910) IBP, IBU, IBV, IBPP, IBF, IUSTVST, IDA 
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WRITE (10, 910) 
WRITE (10, 900) 
WRITE (10, 900) 
WRITE (10, 905) 
WRITE (10, 905) 
WRITE (10, 905) 
WRITE (10,905) 
WRITE (10,905) 
WRITE (10, 905) 
WRITE (10,905) 
WRITE (10, 900) 
WRITE (10, 900) 
WRITE (10, 900) 


ITSSMOO, ISBC, JSK, KSOL, ISKPP 

NX NY , NLE , RX , RY , DXMI N , DYMIN , DXMAX , DYMAX 

SLl, IREADOLD , NSCH , ERM, ERR , ERU, ERV, ERPP , ERF 

PR, PDAM, UDAM, VDAM 

FDAM(l) , FDAM ( 2 ) , FDAM ( 3 ) ,ASW,ES 

TAMB , T S T AR , D AD , Q , E 

FO , XOAMB , VE 0 , VF , EP S 
RSTAR, RS, CS, CP , SIG 
TO, TL,XL, TAU, ASTAR 
QO , QN , QOL , QNL , GLEVEL 

ibeg,jtop,iright,rdxp, DXPMI N , DXPMAX 
ISKIP, I SR, IPARN, VFL1 , VFU1 
MMAX.KPAR, IBC 


•CHECK OUT OVERALL MASS BALANCE IN ELLIPTIC REGION 


FLXBOT=0 . 

FLXTOP=0 . 

FLXRGT=0 . 

FLXLFT=0 . 

RINFASTAR=RSTAR*TSTAR* ASTAR 

E£fc£iS*<..0> *DX (I) /Til. Ojl> 

240 FLXTOP=FLXTOP+V (I , NY) *DX (I ) /T (I, NY, 1) *RINFASTAR 

FLXLFT=FLXLFT+U ( 0 , J) *DY ( J) /T (0, J, 1) *RINFASTAR 
250 FLXRGT=FLXRGT+U (NX, J) *DY (J) /T (NX, J, 1) *RINFASTAR 

FLXNET=FLXBOT+FLXLFT— FLXTOP-FLXRGT 

WRITE (72,*) 'MASS FLUX' 

WRITE (7 2*, 905) FLXBOT , FLXTOP , FLXRGT, FLXLFT, FLXNET 

FORMAT (3(1X,I3),6(1X,F8.4) ) 

FORMAT (IP, 5 (IX, E12. 5) ) 

FORMAT (10 (IX, 13) ) 

FORMAT (3 (IX, F10 . 5, IX, 13, IX, 13) ) , TFFT , 9X 'NET') 

FORMAT (4X,' BOTTOM', 8X, ' TOP', 9X, ' RIGHT', 9X, LEFT , 9X, NET > 

STOP 

SUBROUTINE GRID ( RX,RY, DXMI N, DYMIN, DXMAX, DYMAX) 

COMMON/Ll/DX (0:93) , DY (0:50) , NX, NY, NLE 
REAL RX,RY,DXMIN, DYMIN, DXMAX, DYMAX 
DY (1) =DYMIN 
DX(NLE)=DXMIN 
DX (NLE+1) =DXMIN 
DO 10 J=2 , NY 

10 DY (J) =CVMGT (DY (J) , DYMAX, DY ( J) . LT . DYMAX) 

DO 20 IB=1 , NLE-1 
I=NLE— IB 

20 DX ( I ) =CVMGT^DX^) , DXMAX, DX ( I ) .LT. DXMAX) 

DO 30 I=NLE+2 , NX 

30 DX (I) =CVMGT (DX (I) , DXMAX, DX (I) . LT. DXMAX) 

DX (0) =0 . 

DY (0) =0 . 

DX (NX+1 ) =0 . 


900 

905 

910 

915 

920 
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DY (NY+1 ) =0 . 

RETURN 

END 

SUBROUTINE READIN (TAMB, YOAMB, TSTAR, IREADOLD, IUSTVST ISR 

IBEG, IRIGHT) ' 

COMMON/LI /DX ( 0 : 93 ) , DY ( 0 : 50 ) , NX, NY, NLE 
COMMON/L2/U (0:93, 0:50) , V (0 : 93, 0 : 50) , T (0 : 93, 0 : 50, 1 : 3) , 
/ R(0:93 ' 0:50) ,P(0:93,0:50), PPR(0 : 93, 0 : 50) 

COMMON/ L3/UHAT (0 : 93, 0:50), DU (0:93, 0:50) ,USTAR(0 : 93, 0 ■ 50) 

UAP (93,50) , UAE (93,50), UAW (93, 50) , ' 

rnMMnH/T . 50 > (93, 50) , CRASH, RINFNON 

COMMON/ L4/VHAT (0: 93, 0 : 50) , DV (0 : 93, 0 : 50) , VSTAR(0 : 93, 0 • 50) 

VAP (93, 50) , VAE {93, 50) , VAW(93, 50) , 

“ , VAN (93, 50) , VAS (93,50),VB(93,50) 

COMMON/L8/DXP (0:500) ,DYP (0 : 100) , UPAR (2, 0 : 100) , VPAR (2, 0 • 100 ) 

TPAR (2, 0:100,3) ,TS (0:570,2) , USPAR (0 • 570 ) VSPAR (0 • 570 ) 
COMMON/ L9/IBC, UBC (0:50) , VBC (0:50) , TBC (0 : 50, 3) , VSPAR (0 .570) 

REAL TAMB, YOAMB, TSTAR 
DO 10 M=l,3 


(T (I, J,M) , J=0,NY+1) 
(U (I, J) , J=0, NY+1 ) 


(P (I, J) , J=1 , NY) 

(PPR (I, J) , J=1,NY) 
(USTAR(I, J) , J=0,NY+1) 
(VSTAR(I, J) , J=0,NY) 


DO 10 I=0,NX+1 
10 READ (11,905) 

DO 20 1=0, NX 
20 READ (11, 905) 

DO 30 1=0 , NX+1 
30 READ (11, 905) (V(I, J) , J=0,NY) 

DO 40 1=1, NX 
40 READ (11, 905) 

DO 50 1=1, NX 
50 READ (11, 905) 

DO 52 1=0, NX 
52 READ (11, 905) 

DO 54 1=0, NX+1 
54 READ (11, 905) 

IF (ISR.EQ.l) THEN 

DO 56 I=NLE+1, IBEG+IRIGHT 
56 READ (11, 905) (TS (I , J) , J=1 , 2) 

ENDIF 

IF (IBC.EQ.l) THEN 

READ (11, 905) (UBC ( J) , J=l, NY) 
READ (11, 905) (VBC ( J) , J=1 , NY) 
READ (11 , 905) (TBC ( J, 1 ) , J=1 , NY) 
READ (11, 905) (TBC ( J, 2 ) , J=1 , NY) 
READ (11,905) (TBC (J, 3) , J=1 , NY) 
ELSE 

DO 58 J=l, NY 
UBC ( J) =0 . 

VBC ( J) =0 . 

TBC ( J, 1 ) =0 . 

TBC ( J, 2 ) =0 . 

58 TBC ( J, 3) =0 . 

ENDIF 

IF (IUSTVST. EQ.l) THEN 
DO 60 1=0, NX+1 
DO 60 J=0,NY+1 

USTAR(I, J)=U(I, J) 

VSTAR (I, J) =V (I, J) 

60 CONTINUE 
ENDIF 
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B.C.'S FOR YO- 


CALL DENS (TSTAR) 

SET 

DO 110 J=0 , NY+1 
110 T(0, J,3)=YOAMB 
DO 120 1=0 , NX+1 
120 T (I, NY+1, 3)=YOAMB 

905 FORMAT (IP, 5 (IX, E12. 5) ) 

RETURN 

END 

SUBROUTINE UHATS (PR, TSTAR) 

COMMON/LI /DX (0:93) , DY (0:50) , NX, NY, NLE 

COMMON/L2/U (0:93, 0:50) , V (0 : 93, 0 : 50 ) , T (0 : 93, 0 : 50, 1. 3) , 

R(0:93, 0:50) ,P(0:93, 0:50) ,PPR(0:93, 0:50) 

COMMON/ L3/UHAT (0:93, 0:50) , DU ( 0 : 93, 0 : 50) ,USTAR(0 : 93, 0 : 50) , 

AP (93, 50) , AE (93, 50) , AW (93, 50) , 

AN (93, 50) , AS (93,50),B(93,50), GRASH, RINFNON 
REAL PR, TSTAR, FE, FW, FN, FS, TP , TE, TW, TN, TS, GP, GE, GW, GN, GS , 

DE, DW, DN, DS, GVE, GVW, GVN, GVS, SPDXDY , 

TNORTH,TSOUTH, VOLUME 

; EVALUATE COEFFICIENTS FOR INTERIOR U-CONTROL VOLUMES 

DO 10 1=1, NX-1 

00 FE- J Min,J> . ( 0(I,J)+0(I+1,J) )/2. * OY(J) 

FW= R(T J) * ( U(I, J)+U(I-1/ j) )/2- * DY ( J) 

FN=( R(l' J+1)*DY(J) + R ( I , J) *DY ( J+l ) ) / (DY ( J) +DY ( J+l ) ) 

I +( rTi+I^J+1) *Dy"(j( 2 + R(I+1, J)*DY(J+1) ) / (DY ( J) +DY (J+l ) ) 

* V(I+1,J) * DX(I+l)/2. 

FS= ( R (I , J— 1 ) *DY ( J) + R(I,J)*DY(J-1) ) / (DY ( J) +DY ( J-l 

+( R(I+1,J) *DY ( J— 1 ) + R(I+1, J-l) *DY ( J) ) / (DY (J) +DY (J-l) ) 

_ TP= (T (I J lI (I *DXU + l) + T(I+1,J,1) *DX (I) ) / (DX ( I ) + DX(I + 1)) 

TP— (T (I, J, 1) DXU it T I+2 ' ' * DX(I+1)) /(DX(I+l)+DX(I+2)) 

T (I, J, 1) *DX (I— 1 ) ) / (DX (I ) + DX(I-l)) 


+ 

+ 

+ 

+ 


TE=(T(I+l,J,l)*DX(I+2) 
TW=(T(I-1,J,1)*DX(I) 

TN= (T (I, J+1,1) *DX ( 1+1 ) 

TS= (T ( I , J— 1 , 1 ) *DX (1 + 1 ) 

GP=GAMMAV ( TP , PR, TSTAR) 
GE=GAMMAV (TE, PR, TSTAR) 
GW=GAMMAV (TW, PR, TSTAR) 
GN=GAMMAV (TN, PR, TSTAR) 
GS=GAMMAV (TS, PR, TSTAR) 
DE=2 . * GP*GE / (DX (1+1 ) *GE 
DW=2 . * GP*GW / (DX (I) *GW 
DN=2 . * GP*GN / (DY (J) *GN 
DS=2 . * GP*GS / (DY (J) *GS 
AE (I, J) =DE*A (ABS (FE/DE) ) 
AW (I , J) =DW*A (ABS (FW/DW) ) 
AN ( I , J) =DN*A (ABS (FN/DN) ) 
AS (I, J) =DS*A (ABS (FS/DS) ) 


T(I+1,J+1,1)*DX(I) )/(DX(I) + DX(I+1)) 

T (1+1, J-l, 1) *DX(I) ) / (DX (I ) + DX ( 1+1 ) ) 


+ DX(I+1)*GP) 

+ DX (I) *GP ) 

+ DY (J+l) *GP) * 
DY (J-l) *GP) * 
AMAX1 (-FE, 0. ) 
AMAXl (FW, 0. ) 
AMAX1 (-FN, 0 . ) 


DY ( J) 

DY ( J) 

(DX ( I ) +DX ( 1+1 ) ) / 2. 
(DX (I ) +DX (1+1 ) ) / 2. 


M AMAXl (FS, 0 . ) 

TNORTH =, ( (T (I, J+l, 1) +T (1 + 1, J+l, 1) ) /2 . * DY ( J) 

+(T(I,J,1) +T ( 1+1 , J, 1 ) )/2. * DY (J+l ) ) / (DY (J) +DY (J+l) ) 

TSOUTH= ( (T (I , J— 1 , 1 ) +T ( 1+1 , J— 1 , 1 ) ) /2 . * DY(J) 

+(T(I,J,1) +T (1+1, J, 1) )/2. * DY(J-l) ) / (DY(J) +DY (J-l) ) 

GVE=GAMMAV(T(I+1, J,l) , PR, TSTAR) 

GVW=GAMMAV(T(I, J,l) , PR, TSTAR) 

GVN=GAMMAV (TNORTH, PR, TSTAR) 
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B (I, J) 
-2 


GVS=GAMMAV (TSOUTH, PR, TSTAR) 


20 


- ( GVE/DX (1+1 ) *U (1+1, J) +GVW/DX (I) *U (1-1 , J) )*DY(J)/3. 
*( GVE* ( V ( 1+1 , J) -V ( 1+1 , J-l ) ) — GVW* (V ( I , J) -V ( I , J-l ) ) )/3 

+ ( GVN* (V (1+1, J) -v (I, J) ) -GVS* (V (1+1 , J-l) -v (I , J-l) ) ) 

+GRASH* ( RINFNON* ( DX (I) +DX (1+1) ) 

~ R(I, J) *DX(I)+R(I+1, J) *DX{I+1) ) )*DY(J)/2. 

SPDXDY - ( -GVE/DX (1+1) -GVW/DX (I) ) * DY(J)/3. 

AP (I, J) =AE (I, J) +AW (I, J) +AN (I, J) +AS (I, J) -SPDXDY 
DU (I, J) =DY ( J) /AP (I, J) 

UHAT (I, J) = (AE (I, J) *U (1 + 1, J) +AW ( I , J) *U (1-1, J) + AN ( I , J ) *U (I J+l ) 

- + AS(I,J)*U(I,J-1)+B(I,J))/AP(I,J) 

10 CONTINUE 

DO 20 J=1 , NY 
DU ( 0 , J) =0 . 

DU (NX, J) =0 . 

UHAT (0, J) =U (0, J) 

UHAT (NX, J) =U (NX, J) 

RETURN 
END 

SUBROUTINE VHATS (PR, TSTAR) 

COMMON/L1/DX (0:93) ,DY (0:50) ,NX,NY,NLE 
COMMON/L2/U (0 : 93, 0:50), V(0:93, 0:50), T(0:93, 0:50, 1:3) , 

’ t . R(0:93, 0:50), P(0:93, 0:50) ,PPR (0 : 93,0:50) 

COMMON/L4 /VHAT (0 : 93, 0 : 50) , DV(0 : 93, 0 : 50) , VSTAR (0 : 93, 0 : 50) 

AP (93, 50) ,AE (93, 50) ,AW(93, 50) , 

AN (93, 50) , AS (93,50) ,B (93, 50) 

REAL PR, TSTAR, FN, FS, FE, FW, TP, TN, TS, TE, TW, GP, GN, GS, GE, GW, 

DN, DS, DE, DW, GVN, GVS , GVE, GVW, SPDXDY, 

TEAST, TWEST, VOLUME 

EVALUATE COEFFICIENTS FOR INTERIOR V-CONTROL VOLUMES 

DO 10 1=1, NX 
DO 10 J=1 , NY-1 

FN= R ( I , J+ 1 ) * ( V(I, J)+V(I, J+l) )/2. * DX ( I ) 

FS- R (I, J) * ( V (I, J) +V(I, J-l) )/2. * DX ( I ) 

FE - ( R(I+1, J) *DX(I) + R (I, J) *DX (1+1) ) / (DX ( I ) +DX ( 1+1 ) ) 

* U(I, J) * DY ( J) 12 . 

+ ( R(I+l,J+l)*DXd) + R(I, J+l) *DX(I+1) ) / (DX (I) +DX (1+1) ) 

* U(I,J+1) * DY ( J+l ) /2 . 

FW=( R (1-1 , J) *DX (I ) + R (I, J) *DX (I— 1 ) ) / (DX (I) +DX (1-1) ) 

* U (1-1, J) * DY ( J) /2 . 

+ ( R (1/J+1)*DX(I-1) + R(I-1, J+l) *DX(I) ) / (DX (I) +DX ( I —1 ) ) 

* U(I-1,J+1) * DY (J+l) /2 . 

+ T(I,J+1,1) *DY ( J) ) / (DY ( J) + DY (J+l ) ) 
T (I, J+2, 1) *DY (J+l) ) / (DY (J+l) +DY (J+2) ) 
T ( I , J, 1) *DY ( J— 1 ) ) / (DY ( J) + DY ( J— 1 ) ) 


+ 

+ 

+ 

+ 


TP=(T(I,J, 1) *DY (J+l) 
TN=(T(I,J+l,l)*DY(J+2) 
TS=(T(I, J-1,1) *DY (J) 
TE=(T(I+1, J,l) *DY (J+l) 

TW= (T (1-1, J, 1) *DY (J+l) 

GP=GAMMAV (TP, PR, TSTAR) 
GN=GAMMAV (TN, PR, TSTAR) 
GS=GAMMAV (TS, PR, TSTAR) 
GE=GAMMAV (TE, PR, TSTAR) 
GW=GAMMAV (TW, PR, TSTAR) 
DN=2. * GP*GN / (DY (J+l ) *GN 
DS=2 . * GP*GS / (DY ( J) *GS 
DE=2 . * GP*GE / (DX (I) *GE 
DW=2. * GP*GW / (DX (I) *GW 
AE (I, J) =DE*A (ABS (FE/DE) ) 
AW (I, J) =DW*A (ABS (FW/DW) ) 


T (1+1, J+l, 1) *DY(J) ) / (DY ( J) + DY ( J+l ) ) 
T (1-1, J+l, 1) *DY (J) ) / (DY ( J) + DY (J+l) ) 


+ DY (J+l) *GP) 

+ DY ( J) *GP ) 

+ DX (1+1 ) *GP) * 
+ DX(I-l) *GP) * 
+ AMAX1 (— FE, 0 . ) 
+ AMAX1 (FW, 0 . ) 


DX ( I ) 

DX(I) 

(DY ( J) +DY ( J+l ) )/2. 
(DY ( J) +DY (J+l) ) /2 . 
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S&SSKSSSSi! : JKKTiV 

- “5^ cx " "7 1 r Ji ’ t/2. * DX(I-I) ) / (DX (i) +dx ( i 1 > ) 
GVN=GAMMAV (T (I, J+l, 1) , PR/ TSTAR) 

GVS=GAMMAV (T (I, J/l) f PR/ TSTAR) 

GVE=GAMMAV (TEAST/PR, TSTAR) 

GVW=GAMHAV(TWEST, PR, TSTAR) . * DX (I)/3. 

V.{ S5K58TiinS’,r-i ;'»• 

SPDXD J «.i, 3SCT.ZS35 £ j) — SPDXDY 

+ AS (I, J) *V (I, J-l) +B (I, J) ) / AP l 1 ' J) 

10 CONTINUE 

DO 20 1=1, NX 
DV (I , 0 ) =0 . 

DV (I , NY) =0 . 

VHAT (I, 0) =V (1,0) 

20 VHAT (I, NY) =V (I, NY) 

RETURN 

SUBROUTINE PRESSURE (ERP, PDAM, IBP) 

COMMON/L1/DX(0:93),DY(0:50),NX,NY,NLE 

COMMON/L2/U (0:93, 0:50) 'p |q| gg' gl jqj ^ ppr"( 0^93 j 0: 50) 

COMMON/ L3 /UHAT ( 0 : 93 , 0 : 50) , DU ( 0 : 93 , 0 : 50 ) ^USTAR (0:93,0.5 ) , 

: San (II' 50) ' UAS (93 f , 50) ', UB (93, 50) , GRASH, RINFNON 

COMMON/ L4 /VHAT (0:93,0:50) , DV (0:93,0:50) , VSTAR(0 : 93, 0 : 50 ) , 

VAP (93, 50) , VAE (93, 50) ,VAW (93, 50) , 

VAN(93, 50) ,VAS (93, 50) ,VB(93, 50) 

COMMON / L5 / AP (93, 50) , AE (93, 50) ,AW (93, 50) , AN(9! 3^50) , 

AS(93,50),SC(93,50) ,SP(93,50),B(93,50) 

COMMON/L6/AA(0:93),BB(0:93),CC(0:93),DD(0:93),XX(0. 

S ^HS, RESP ERP, PDAM, 

c : a? SiS« cm. VOLUHES 

DO 10 1=1, NX 

DO 10 J=1 , NY , 4 x 

FFE=DX ( 1 + 1 ) / (DX ( I ) +DX ( 1+1 ) ) 

FFW=DX (I— 1 ) / (DX(I)+DX(I-1) ) 

FFN=DY (J+l) / (DY ( J) +DY (J+l) ) 

FFS=DY (J-l ) / (DY ( J) +DY (J-l) > 

RHE= (R (I , J) *FFE+R (1+1, J) * (l.-FFE) ) 

RHW= R I J)*FFW+R(I-1,J)*(1.-FPW)) 

Sn= (R (I , J) *FFN+R (I , J+l ) * (1 • -FFN) ) 

RHS= (R (I , J) *FFS+R (I, J-l) * (l.-FFS) ) 

AE (I , J) =RHE*DU (I , J) *DY ( J) 

AW (I, J)=RHW*DU (1-1, J) *DY (J) 

AN (I , J) =RHN*DV (I , J) *DX (I ) 

AS (I, J) =RHS*DV(I, J-l) *DX ( I ) 

B (I, J) = (RHW*UHAT (1-1, J) — RHE*UHAT (I, J) ) DY ( J) 
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10 


20 


30 

40 


f (RHS*VflAT (I, J-l) — RHN*VHAT (I, J) ) *DX (I) 

CONTINUE ' ~ AE ( 1 ' J> +AW ( 1 ' J> +AN ( 1 ' J) +AS ( 1 ' J) 

CALL OLDS (POF , P, NX, NY) 

DO 80 IBAILOUT=l, IBP 
CALL OLDS (PO,P,NX,NY) 

~DO~ 40 ~I=l"ra~" SWEEP USING VERTICAL LINES— 

DO 20 J=1 , NY 

AA(J)=AP(I, J)/PDAM 
BB ( J ) =AN ( I , J) 

CC ( J) =AS (I, J) 

DD(J)=AE(I, J)*P(I+l,J)+AW( I ,J)*p (I - 1 , J)+ B (i j) 

+ ( 1 . -PDAM) /PDAM* AP ( I , J) *POF ( I , J) 

IF (I.EQ.1|_THEN ' ' 

cc"(ny)"=o ET pressure ar bitrarily=o AT TOP LEFT POINT 

AA (NY) =1 . 

DD (NY) =0 . 

ENDIF 

CALL TDMA ( 1 , NY ) 

DO 30 J=1 , NY 
P (I, J) =XX (J) 

CONTINUE 


50 


60 

70 


80 

90 


nn 7 n t i SWEEP USING HORIZONTAL LINES 

/ U J— 1 f NY 

DO 50 1=1, NX 

AA (I) =AP (I, J) /PDAM 
BB (I) =AE (I, J) 

CC (I) =AW (I, J) 

DD (I) * P <1/ J+l) +AS (I, J) *p (I, J— 1) + B (I, J) 

+ (1 • -PDAM) /PDAM*AP (I, J) *POF (I, J) 

IF (J.EQ.NTO THEN ' 

a ^“~““~ SET PRESSURE ARBITRARILY=0 AT TOP LEFT POINT ( 

BB (1) =0 
DD (1 ) =0 . 

ENDIF 

CALL TDMA (1, NX) 

DO 60 1=1, NX 
P (I, J)=XX(I) 

CONTINUE 

IF (ABS^SPKLT^Sj^O^^SO^ 3501 ^ '' IDIF ' JDIF > 

CONTINUE ~ ' v 

CONTINUE 

RETURN 

END 

SUBROUTINE USTARS (ERU, UDAM, IBU, PR, TSTAR) 

COMMON/ LI /DX (0 : 93) ,DY (0 : 50) , NX, NY, NLE 

COMMON/L2/U(0:93, 0:50), V(0:93, 0:50), T(0-93 0-50 1-3) 
R(0:93,0:50),P(0:93 0:50 PPR(0-9?n ' 

COMMON/L3/UHAT ( 0 : 93, 0 : 50) , DU (0:93,0:50), USTAR ( 0 • 93 0‘50 
AP (93,50) , AE (93, 50 ) , AW (93,50), ' * 5 ° 

CQMMON/L6/J^0?^ 

COMMON/L9/IBC,UBC(0:50),VBC(0:50),TBC^-56 3 ’ 

real USTAROF (0:93,0:50), USTARO (0:93,0:50), RESU, ERU, UDAM 


DOMAIN 


DOMAIN 



o o 
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C 


CALL OLDS (USTAROF , USTAR, NX, NY) 

DO 80 ibailout=i,ibu 

CM.L o MS1 USTMO,US I lise5 - 


DO 40 1=1, NX-1 
DO 20 J=1 , NY 

AA ( J) =AP (I , J) /UDAM 
BB(J)=AN(I, J) 


20 


+ (l^-UDAM) /UDAM*AP (I, J) *USTAROF (I, J) 

— + U ' 


■IMPOSE BOUNDARY CONDITIONS 


- bottom--^ (NLE>1)) then 

AA ( 0 ) =1 • 

BB (0) =1 ■ 

DD ( 0 ) =0 . 

ELSE 
AA ( 0 ) =1 • 

BB (0 ) =0 . 

DD { 0 ) =U (1,0) 

ENDIF 
C — TOP — 

CC (NY+1) =0 . 

AA (NY+1 ) =1 • 

DD (NY+1) =1 • 

CALL TDMA (0 , NY+1) 

DO 30 J=0 , NY+1 
30 USTAR (I , J) =XX ( J) 

4 0 ^CONTINUE gwEEp H0RIZ0NT AL LINES 

DO 70 J=1 , NY 
DO 50 1=1, NX-1 

AA(I)=AP(I,J)/U DAM 
BB (I) =AE (I, J) 

£ a! a : *UST« (I, r n ♦« ,i a, -ustar u . j-u +b (i f 

I ‘ a !Si2o ToSSfii' <U> *usT»or ,1.,. 

-—impose boundary conditions — 


50 


c — LEFT — 

AA (0) =1 • 

BB (0) =0 . 

DD (0) =1 • 

C — RIGHT — 

CC (NX) =1 - 

OT (S) =DX (NX) *UBC ( J) /GAMMAV (T (NX, J, D , TSTAR) 

CALL TDMA (0, NX) 

DO 60 1=0, NX 
60 USTAR (I, J) = XX (I) 

70 SE'EW (USTAR, USTARO, NX HY , BESU, - U-STAR •.IDIF.JPIF, 
IF (ABS (RESU) .LT.ERU) GO TO 90 
80 CONTINUE 
90 CONTINUE 
RETURN 
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20 


END 

SUBROUTINE VSTARS (ERV, VDAM, IBV PR TSTflRi 
COMMON/ LI /DX <0 : 93) ,DY <0 • 50) Sr' It l ' 

CALL OLDS (VSTAROF; VSTAR'ra T w? (0 : 93 ' ° : 50) ' ER V, VDAM 
DO 80 IBAILOUT=l , IBV K,NX ' NY) 

CALL OLDS (VS TARO, VSTAR, NX, NY) 

DO 4 0 I=i7nx” SWEEP USING VERTICAL LINES 

DO 20 J=l, NY-1 

AA ( J) =AP (I, J) /VDAM 
BB ( J) =AN (I, J) 

CC (J) =AS (I, J) 

DD ,J) - “t n ♦»«.» 



IF (I.LE.NLE) THEN 
AA ( 0 ) =1 . 

BB (0) =0 . 

DD (0) =0. 

ELSE 

AA (0) =1 . 

BB (0) =0 . 

DD ( 0 ) =V (1,0) 

ENDIF 
C — TOP — 

CC (NY) =1 . 

AA (NY) =1 . 

DD (NY) =0. 

CALL TDMA (0, NY) 

DO 30 J=0, NY 
VSTAR (I, J) =XX ( J) 

CONTINUE 

Do"o"j:r;^:r SWEEP USING horiz ° nt al lines 

DO 50 1=1, NX 

AA ( I ) =AP ( I , J ) /VDAM 
BB (I ) =AE (I, J) 

CC (I) =AW (I, J) 

DD zngxn r «.«»««!.« 




30 

40 


50 


C — LEFT — 

AA(0)=1. 

BB (0) =0 . 

DD ( 0 ) =0 . 

C — RIGHT — ‘ 

CC (NX+1 ) =1 
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AA (NX+1) =1 . 

TMEAN=(T(NX, J,1)+T(NX+1, J,l) ) /+. 

DD (NX+1) =DX (NX) /2 . *VBC ( J) /GAMMAV (TMEAN, PR, TSTAR) 

CALL TDMA ( 0 , NX+1 ) 

DO 60 1=0, NX+1 
60 VSTAR(I, J) =XX (I) 

70 CAL^FINDDIF (VSTAR, VSTARO, NX, NY, RESV, ' V-STAR ' , IDIF, JDIF) 

IF (ABS (RESV) .LT.ERV) GO TO 90 
80 CONTINUE 
90 CONTINUE 
RETURN 
END 

SUBROUTINE PPRIME (ERPP , IBPP , ISKPP ) 

COMMON/L1/DX(0:93),DY(0:50),NX,NY,NLE 

COMMON/L2/U(0:93,0:50),V(0:93,0:50,T (0:93, 0.50, 1.3), 

r(0*93 0 * 50) , P (0 : 93/ 0 : 50 ) , PPR (0 : 93, 0 . 50 ) 

COMMON /L3/UHAT (0:93,0:50) , DU ( 0 : 93, 0 : 50) ,USTAR (0 : 93, 0 : 50 ) , 

TJAP (93.50) , UAE (93, 50) , UAW (93, 50) , 

UAN (93,50) , UAS (93, 50) ,UB(93, 50 ) , GRASH, RINFNON 
COMMON/L4/VHAT (0:93, 0:50) , DV ( 0 : 93, 0 : 50) , VSTAR (0 : 93, 0 : 50 ) , 

VAP (93,50) , VAE (93,50) ,VAW(93,50) , 

VAN (93,50) , VAS (93, 50) ,VB{93,50) 

COMMON / L 5 / AP (93,50) ,AE (93, 50) , AW (93, 50) ,AN(93, 50) , 

AS (93, 50) , SC(93, 50) , SP (93, 50) ,B(93, 50) 

COMMON/ L6/AA (0 : 93) ,BB(0:93) , CC ( 0 : 93) , rhf^RHW^ RHN RHS 
REAL PPRO (0 : 93, 0 : 50 ) , FFE, FFW, S ^> P 

- FE FW, FN, FS, DE, DW, DN, DS, GE, GW, GN, GS, GP, RESPP, ERP 

EVALUATE COEFFICIENTS FOR INTERIOR CONTROL VOLUMES 

DO 10 1=1, NX 
DO 10 J=1 , NY 

FFE=DX ( 1+1 ) / (DX ( I ) +DX ( 1+1 ) ) 

FFW=DX ( I— 1 ) / (DX (I) +DX (1-1) ) 

FFN=DY (J+l) / (DY (J) +DY (J+l) ) 

FFS=DY { J — 1 ) / (DY (J) +DY (J-l) ) 

RHE= (R (I , J) *FFE+R(I+1, J) * (l.-FFE) ) 

RHW= (R (I , J) *FFW+R ( I— 1 , J) * (l.-FFW) ) 

RHN=(R(I, J)*FFN+R(I, J+l) * (l.-FFN) ) 

RHS= (R (I , J) +FFS+R (I , J - l) * (l.-FFS) ) 

AE (I , J) =RHE*DU (I, J) *DY ( J) 

AW (I, J)=RHW*DU (1-1, J) *DY (J) 

AN (I, J) =RHN*DV (I, J) *DX ( I ) 

AS (I, J)=RHS*DV(I, J-l) *DX(I) 

B (I , J) = (RHW*USTAR (1-1 , J) -RHE*USTAR (I, J) ) *DY ( J) 

+ (RHS*VSTAR (I, J-l) -RHN*VSTAR(I, J) ) DX(I) 

AP ( I , J) =AE ( I , J) +AW ( I , J) +AN ( I , J) +AS ( I , J) 

10 CONTINUE 

DO 80 IBAILOUT=l, IBPP 

CALL OLDS (PPRO, PPR, NX, NY) 

SWEEP USING VERTICAL LINES 

DO 40 1=1, NX 
DO 20 J=1 , NY 
AA(J)=AP(I, J) 

BB ( J) =AN (I, J) 

DD ( J) =AE ( I J) *PPR (1 + 1, J) +AW (I, J) *PPR(I-1# J) +B (I, J) 

TOP RIGHT POINT OF DOMAIN 


20 
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50 


30 


40 


50 


CC (NY) =0 . 

AA(NY)=1. 

dd(ny)=o. 

ENDIF 

CALL TDMA (1,NY) 

DO 30 J=l, NY 
30 PPR (I , J) =XX ( J) 

40 CONTINUE 

SWEEP USING HORIZONTAL LINES 

DO 70 J=1 , NY 
DO 50 1=1, NX 
AA(I)=AP (I, J) 

BB (I) =AE (I, J) 

CC (I) =AW (I, J) 

DD (I) =AN (I, J) *PPR (I, J+l ) +AS (I, J) *PPR(I, J-1)+B(I J) 

IF ((J-EQ.NYKAND. (ISKPP.EQ.O)) THEN U+B(I,J) 

CC(NX)”=0 ET p_prime arbitrarily=o at top right point of domain 

AA (NX) =1 . 

DD (NX) =0 . 

ENDIF 

CALL TDMA (1, NX) 

DO 60 1=1, NX 
60 PPR (I , J) =XX (I ) 

70 CONTINUE 

CALL FINDDIF (PPR, PPRO, NX, NY, RESPP, ' P-PRIME ' IDIF TTVTF1 
IF (ABS (RESPP KLT.ERPP) GO TO 90 ,IDIF,JDIF) 

80 CONTINUE 
90 RETURN 
END 

SUBROUTINE CORVEL (PR, TSTAR) 

COMMON/ LI /DX (0:93) , DY (0:50) , NX, NY, NLE 
COMMON/L2/U (0:93, 0:50), V(0: 93,0: 50), T{0: 93, 0:50, 1:3) 

R (0 : 93, 0 : 50) , P (0 : 93, 0 : 50) , PPR(0 • 93 0-50) ' 

^COMMON/L3/UHAT ( 0 : 93, 0 : 50 ) , DU (0:93,0:50) , USTAR (0:93,0:50) 

UAP (93, 50) ,UAE (93, 50) ,UAW(93, 50) , 

50) , UAS (93, 50) , UB (93, 50) GRASH r tnfmom 
COMMON/L 4/VHAT < 0 : 93, 0 : 50 ) ,DV ( 0 ! 93, oTiof, hTSS) 

VAP (93, 50) , VAE (93, 50 ) , VAW (93,50) , ' ' ’ ' 

VAN (93,50) , VAS (93,50), VB (93,50) 

DO^O^-l^ 80 ' 1180 (0:50) ' VBC <°:50) ,TBC(0:50,3) 

DO 10 J=1,NY 

nn J)+DU(I ' J > * (PPR(I, J)-PPR(I+1, J) ) 

Z U 1 = 1, NX 

DO 20 J=1 , NY— 1 

V(I, J) =VSTAR(I, J) +DV (I, J) * (PPR (I, J) -PPR (I, J+l ) ) 

ADJUST BOUNDARY VELOCITIES------ 

DO 30 J=1 , NY 

TMEAN= (T (NX, J, 1 ) +T (NX+1 , J, 1 ) ) /2 
GV1=GAMMAV (T (NX, J, 1 ) , PR, TSTAR) 

GV 2 =GAMMAV (TMEAN, PR, TSTAR) 

U (NX, J) = (U (NX, J) +U (NX-1 , J) +DX (NX) *UBC { J) /GV1 ) /2 

D0^4 0^1=1 ^NX <V (NX+1 ' J) +V (NX ' J) +DX (NX > / 2 • *VBC ( J) /GV2 ) /2 . 

V ( I , NY ) =V ( I , NY— 1 ) 

DO 50 1=0, (NLE-1) 

U(I,0)=U(I,1) 


10 


20 
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RETURN 

END ~ trn cpp FDAMr M* TSTARf lor } 

mS / m"o 5 93K d? JO /so! (0:93,0:50,1:3), 

^OHMOH/l 2 /mO = 33,0 ;| 0) ^ S|;0:50,, fo ,;-».0:50, 

_C0MM0N/L5/AR (93, 50, ; £ ■ ’ Si > (|3 50, ,B,|3, 50, ^ 

Co m OS/16/«(0:93),BB(0:9^,CC(0 93K (0;50 3) 

SISKS 35 :*“^ hi ^1°, •££%&.»»»■ 

_ TAMB, DA, Q-E, ERF, FDAM (3) , PEF, FO^ TO NOFF , TE1 , 

tonoff=o . 

IF (M.EQ.l) THEN 
PRESP=0 . 

PRESC=Q 

tonoff=i . 

ENDIF 

IF (H.EQ.2) THEN 
PRESP=“1 • 

PRESC=0 . 

ISPEC=3 
LE(M)=1- 
YI (M) =1 • 

ENDIF 

IF (M.EQ.3) THEN 
presp=-fo 

PRESC=0 . 

ISPEC=2 
LE (M) =1 • 

YI (M) =0 . 

.-..^evaluate coefficients for interjor control volume 
DO 10 1=1, NX 

D0 FFE=DX^I+1 > / (DX (X)+DX(I+1) > 

f ™=dx (i-D / +DX 1 H! 

FFN=DY (J+l) / (DX (J) +DX ( J+l) > 

fS=DY(J-1)/(DX(J)^Y(J-1)) 

FE=(R(I,J>* FFE+R(I+ ^ *„ 

FW= (R(I, J) * FFW+R( ,t"t+i! * (1 
FN=(R(I, J) *FFN+R(I, J+l | ( ’_ffS) ) 

fs= (r (i» j) ’ 

GP=GAMMAT (T (I, j, 1) -M, I s *™ ) 

GE=GAMMAT (T <1+1 , J' 1) ' JI^AR) 

gw=gammat (T (1-1/ J» 1 »*J' I s ™ R 

GN=GAMMAT (T (I , j+ 1, 1 ) ' M, 

GS=GAMMAT *GE +’ DX (1+1) *GP) 

sps- : ss / KimS + 

SS-2 * GP*GN / (DX(J)*GN 

Sc=2 * GP*GS / (DY(J)*GS 
AE(I j) =DE*A (ABS (FE/DE) ) + 

AW(l',J)=DW*A(ABS(FW/DW)) + 

an (I , J) =DN*A (A BS (FN/DN) ) 

AS(I,J)=DS*A(ABS(FS/DS) ) + 


. — FFE) ) 
,-FFW) ) 
. -FFN) ) 


U(I, J) 
uu-i, J> 

V(I, J) 
V(I, J-i) 


* DY (J) 

* DY ( J) 

* DX ( I ) 

* DX ( I ) 


+ DY(J+1)*GP) 
+ DY (J-l) *GP ) 
AMAX1 (-FE, 0 . ) 
AMAX1(FW,0.) 
AMAX1 (-FN, 0 . ) 
AMAX1 (FS, 0 . ) 


DY ( J) 
DY ( J) 
DX (I) 
DX ( I ) 



o o 


134 


B (I, J) =SCC*DX (I) *D^'(J) 2 U,J ' 3) * ARRHEN 

1 0 CONTINUE ' J) ~ AE ( 1 ' J> +AW ( 1 ' J) +AN ( 1 ' J) +A S ( X, J) -SPP *DX ( I ) *DY ( j) 
CALL OLDS3 (TOLF, T, NX NY mi 
DO 80 IBAILOUT^Bf' ' M) 
c CALL^OLDS3 (TOL, T,NX,NY,M) 

DO~4 0 ~l=l ~nx SW EEP USING VERTICAL LINES 

DO 20 J=1 , NY 

AA ( J) =ap (I , J) /FDAM (M) 

BB ( J) =AN (I, j) 

?n <-C(J)=AS(I,J) 

— BOTTOM IMPOSE BOUNDARY CONDITIONS l_l_i 

IF (I.LE.NLE) THEN 
AA(0)=1. 

BB (0) =1 . 

DD (0) =0 . 

ELSE 

• x , , M . TSTAR) , 

BB ( 0 ) =1 . -TONOFF 

eJ£i? )=TE1 * YI (M)+T <Ia0,1)*tONOFF 
C — TOP — 

CC (NY+1 ) =0 . 

AA (NY+1) =1 . 

r'M? ( f +1)=T(I ' NY+1 ' M) 

CALL TDMA (0, NY+1 ) 

DO 30 J=0, NY+1 

30 T(l, J,M)=XX(J) 

40 CONTINUE 

DO 70 J=1~NY~~ SWEEP USING HORIZONTAL LINES 

DO 50 1=1, NX 

AA (I) =AP (I, J) /FDAM (M) 

BB ( I ) =AE ( I , J) 

CC (I ) =AW (I, J) 

C— — LEFT— — 1OT0S E BOUTOARY c6t^T™s---~I:f:”l 

AA (0) =1 . 

BB (0) =0 . 

C- R IGHT-- (0, " r(0 ' J ' M) 

CC (NX+1 ) =1 . 

AA (NX+1 ) =1 . 

TMEAN= (T (NX, J, 1 ) +T (NX+1 J l)\/o 
«£«’,;££!' ’^THEAN.M.TSTAR, 
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DO 60 1=0 , NX+1 
60 T(I» J,M)=XX(I) 

10 SE^FS ,*.^"-?"-'" 1 '•«> 

IF (ABS (REST) -LT. ERF) GO TO 90 
80 CONTINUE 
90 CONTINUE 
RETURN 

ComON/LmX{0:93l,DF<°^OWm^l«jtILE 3 o. 50/ i ; 3) , 

COMMON/L2/U (0: 93 ; 0 : |0) 'p jgj* 0 I 50 , ^PPR(0: 93, 0:50) 

REAL tstar 
DO 10 1=0, NX+1 
DO 10 J=0,NY+1 
10 R (I , J) =TSTAR/T ( I , J, 1 ) 

return 

SS”"93;:bb«0 : 93, ,CC <0 : 93, ,BO <0 : 93, ,XX <0 : 93, 

C-THE^ll OF T»’e U <11 XX (I > -BB <1, XX (I* I ♦= (I , « > + ™ 111 

c— FORWARD SUBSTITUTION— 

CC(N1)=0. 

BB (N2) =0 . 

xo Sc -occ. — «-« . 

C— BACK SUBSTITUTION— 

XX (N2) =QQ (N2) 

DO 20 II=N1 , N2-1 
t =N2-1+N1““H 

20 XX I (I)=PP (I)*XX(I+D+QQ (I) 
return 

SUBROUT I NE WRITE1 (X, TX, NXO ,NXT, NFORMAT, 

REAL X (0 : 570) 

IF (TX EQ . 0 . ) TX— 1 . 

IF „RS I 9io?'V ( xamx,. I .Hxo. a xT, 

WRITE (13, 905, < (X (I, /TX) , I-NX0.NXT, 

endif 

C 900 FORMAT (5 (IX, F7. 4) ) 

905 FORMAT (IP, 5 (IX, E12 . 5) ) 

return 

end 

C SUBROUTINE WRITE2 (X, NXO , NXT, NY 0, NTT, NFORMAT, 

REAL X(0:93 f 0:50) 

IF (NFORMAT . EQ . 1 J THEN 
DO 10 I=NX0 , NXT t-myO NYT) 

10 WRITE (13, 900) 

ELSE 

DO 20 1= NXO , NXT Krvr, myt^ 

20 ^ WRITE (13,905) (X (I, J> . J-*«0.NYT> 
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900 

905 


endif 

FORMAT (5 (IX, F7. 4)) 
FORMAT (IP, 5 (lx, E12. 5 )) 

RETURN ' 

END 


10 


20 


900 

905 


10 

20 


REAL X (0:93, 0-50 3) '' NX0,NXT ' NY0 ' NYT /NFORMAT) 
DO 5 M=2, 3 ’ ' 

DO 5 I=NX0, NXT 
DO 5 J=NY0, NYT 

i CONTrHUE' J,M) ■ LT - 1 - e-7 8) X(I,j,h). 0 . 

^DO ^ ^ 

DO 10 I=NX0, NXT 

ELSE ITE <13 ' 9 ° 0) <X (I ' J ' M > ' J=NY0, NXT) 

DO 20 M=1 , 3 
DO 20 I=NX0, NXT 

ENDIF 1 ^ (13 ' 905> (X (I ' J ' ' J=NY0, NYT) 

FORMAT (5 (IX, F7. 4) ) 

FORMAT (ip 5 (lx £12 5)) 

RETURN " 

END 

REAL X^of^O^^ 4 (X ' '^ X ®' NXT, NY0, NYT, NFORMAT) 
!F (NFORMAT.EQ.i) THEN 
DO 10 I=NX0, NXT 

else ITE(13,900) (X(I ' J )/J=NY0,NYT) 

DO 20 T=NX0, NXT 

endif ITE<13, 905) (X(I ' » j= nyo,nyt) 


900 

905 


10 


FORMAT (5 ( IX, F7. 4) ) 

FORMAT (IP, 5 (lx £12 5)) 

RETURN 

END 

Sa^aSHc”^^ VARNAME, IDIF JD1F) 
CHARACTER VARnWiO ' ARR2(0:93 ' 0 * 50 > .MAXDIF, LOCDIF ' 1F) 

maxdif=o . 

idif=o 

jdif=o 

do 10 1=0, NX+1 
DO 10 J=0, NY+1 

LOCDIF=ARRl (I, J) -ARR? (T t \ 

IF Mi5DIF=LSl^ GT - ABS( ^ DIF >> THE ” 

IDIF=I 

JDIF=J 

ENDIF 

CONTINUE 
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. 13) 


WRITE (12,900) '' I3) 

900 FORMAT (' DELTA , A, ' 

RETURN 

CHARACTER VARNAME*10 
MAXDIF=0 . 

IDIF=0 
JDIF=0 

' DO 10 1=0, NX+1 

maxdif=locdif 
IDIF=I 
JDIF=J 
ENDIF 

10 SScaI.900, 

900 FORMAT (' DELTA ,A,I2, 

RETURN 

|S°S£® 93°0 = 50, 

DO 10 1=0 ,NX+1 
DO 10 J=0 , NY+1 

10 OLD (I, J) “CURRENT (I, J) 

RETURN 

END 

° SUBROUTINE OLDS3 (OLD, CURRENT, NX, NY,M) 

S°oS(0: 93,0:50,3), CURRENT (0:93,0. 50,3) 

DO 10 1=0 , NX+1 

10 D ° 0 LD “CURRENT (I, J,M) 

return 

DO 10 1=0, NX+1 

10 . + — < ,. 0 , 

return 

END 

c srS"“ 

10 "S» «««»».>.» * (i - fract> * oi - dval a • M) 

return 

function CVMGT (AAA,BBB,CCC) 
real aaa,bbb 
logical ccc 
if (CCC) THEN 
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CVMGT=AAA 

ELSE 

CVMGT=BBB 

ENDIF 

RETURN 

END 

FUNCTION A (P) 

REAL~X~P PECLET " FUNCTI0n PREFERRED 
X= (1 . — 0 . l*p ) **5 
IF (X.GT.O.) THEN 
A=X 
ELSE 
A=0. 

ENDIF 

RETURN 

END 


IN SIMPLE ALGORITHM- 


C FUNCTION^GAMMAT (XXX, M, TSTAR) 

COra >“CTlVIT r , GAMMAT, at T-XXX 

GAMMAT=XXX / TS TAR 

RETURN 

END 

C 


GAMMAV (*XX, PR, TSTAR) 

real xxx,prJstar ATES VISC0SITY ' GAMMAV, at T=XXX AND GIVEN PR— 

GAMMA v=p R * XXX / TS TAR 

RETURN 

END 

^SUBROUTINE SOLID (EPS ,VFL, VEU. VE. VEO, UBUOT. 

S :|SJ ; i 5 J • 3 . - 

^COMMON / LS/DXP^(0^500) ,DYP^((h Joo?|upar( 2^0: 100) ^EPAR^Jf 0^100) 

w ' ™ - xndanp. flsw , UREF 

RR=rs*cs/rstar/cp 

f=f^* EPS *T0** 3 /RSTAR/CP 

C^CP/2./CS CS) * TL/T0 ~ XL/CS/T0 

E=ES 

VFLINI=VFL 
VFU INI =VFU 
THKNEW=1 . 

THKOLD=0 . 

To'eeriiiir;?^™ mmm routine 

THKOLD=THIUIEW KOLD/THKNEW) ' LT ' 0,01) GOTO 999 

VFL=VFLINI 

VFU=VFUINI 
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10 


20 


MM=0 

DO 10 IK=NLE+1,NX2+1 
TOV(IK)=TS (IK,1) 

<TS <1. 2, *TS <I,D . MTS (1 . 2, -TS (1. 1> ) /« «> /TS " R 

~ ,!?.SSr i ;!Su.i> » * (« ti. ^ /TSTOR 

CONTINUE 
30 VF= (VFL+VFU) /2 . 

MM=MM+1 

IF (MM.GT.MMAX) THEN _ , VF 

WRITE (*,*) 'MM EXCEEDED, VF - /VF 

ZZZ=1. 

GO TO 999 
ENDIF 

UREF=VEO+UBUOY-VF 

XR=AS TAR/ UREF 
ri=rr*vf/uref 
S1=S/UREF 

THICK (NX2+1 ) =TAU/XR 
DO 40 J=NLE+1 , NX2 

Sx=S/VF*( / ^XNDAMP+XODAMP ) 

TF (I LT IBEG) DDX-(DX(I)+DX(I+l))/2. 

J !F I'.GE.IBEGi DDX-DXP (I-IBEG+l) 

THICK ( I ) =THICK (1+1 ) -DHDX*DDX 
IF (THICK(I).LT.O.) THEN 
VFL=VF 
GO TO 30 

/GDX-OHD X MBN*C.TS(I + l, 1> ) . 

riUS«(n/DDxl0HDX.CI + 2..Sl.TS ( I + l. l )**3 
TS(I 1) = (QYA+AA+BB) /CC 

" ‘(isd'll.W.l-l TSII.11-1. 

CONTINUE 


40 


50 


888 


TINUEi 

IF (ABS (THICK (NLE+1) ) . LT . 

IF (THICK (NLE+1) .GT . 0 .DO) 

VFU=VF 
GO TO 30 
ENDIF 
CONTINUE 

tx=tau/xr 

THKNEW=THICK (NLE+1) 

IF (KSOL.NE.O) GO TO 999 


(1 ,D-4*TAU/XR) ) 
THEN 


GO TO 50 


900 FORMAT (IX, OP/ 13/ IP/ 4 E g 2) 

99 5 9 SSSaKoSr^vr: (TH IC K (NLE+1) /TX) 

CONTINUE 

RETURN 

END trEG JTOP , IRIGHT , RDXP , DXPMIN , DXPMAX) 

93) ,0 : 93, 0 : 50, 1 : 3) , 

^COMMON/L2/U (0:93, 0:50) 'p jgj* g! 50 ) (pPR(0 : 93,0:50) 
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10 

20 

30 

33 


35 

40 

50 


10 


20 

25 

C**** 

c**** 

c**** 

30 


-COMMON/L8/DXP^0^500) ,DYP (0^100) ,UPAR(2, 0:100) ,VPAR(2, 0 : 100) , 

IF (MF.BQ.l) GO TO 35 °' 3) ' TS «> • 570, 2) , USPAR <0 : 570) , VSPAR <0 • 570 ) 

FFC=1 . 

DO 10 J=0,NY 
DYP(J) =DY ( J) 

DO 20 J=NY+1 , NY+1+ j T0P 
DYP ( J) =DY (NY) 

DO 30 I=IBEG, NX 

DXP (I-IBEG+1) =DX (I) /2 . +DX (I+ll /? 

dxp(nx +2 -ibeg)=dxpmin DX(i+1)/2 - 
3 i=nx+3-ibeg,iright 

DXP (I) -DXP (1-1) *RDXP 

--^l^^^P^J'DXPMAX^XPfl, .LT DXPMAY> 

IF ( IBEG . ^ S+l ) V FFC-0 F0R PARABOLIC 'cALCUllTION 

DO 40 J=0, NY+1 

TPAR ( 1 ' J 1 )=T aBEcf jfl f } +U ( IBEG_1 ' J > > / d.+FFC) 

TPAR (1, J, 2) =T (IBEG, J, 2 ) 

TPAR(1, j, 3) =T (IBEG, J, 3) 

VP AR ( 1 , 0) =V (IBEG, 0) 

DO 50 J=1,NY 

^S;S«7SS!iS!S.' I “' i J ' 1 ’ ’ n 

DO 60 J=NY+2,NY+l+jTOP 
UPAR (1 , J) =UPAR (1 , NY+1 ) 

VPAR ( 1 , J) =VPAR ( 1 , NY+1 ) 

TPAR(1, J, 1) =TPAR(1,NY+1 1) 

TPARU, J, 2 )=TPAR(1,NY+1 2 

TPAR (1, J , 3) =TPAR ( 1 , NY+1 ^ 3 ) 

CONTINUE ^ 

RETURN 

END 

-SUBROUTINE SETSUR (IBEG, JTOP, I RIGHT, ES, ASW, RS 
COMMON/L1/DX (0 : 93) ^Y^O • 50 f ^N^my^ 1 *' ISR) 

j : r|«^”65;S?5Si«' 

» «»..- a T;r<g*& i sS' 3> - ” i T «, 

DO 10 I=NLE+1, IBEG 
TS(I,1)=T(I,0,1) 

TS (1,2) =T (1,1, 1) 

DO 20 I=IBEG+1, IBEG+IRIGHT 

TS (1, 1) =TS (IBEG, 1) -(TS (IBFr ii \*,- r 

TS (1,2) =TS (IBEG, 2) - (TS (IBEG , 24 ~i * l * j*~ IBEG ) / d • *IRIGHT) 

DO 30 I=IBEG, IBEG+IRIGHT ' _1 ’ d-IBEG) / (1 . BRIGHT) 

•*.*”**ss ”r™ssss*:s**************«*****..... 

******************** Mo-rsiwx 


60 


VSPAR(I) =EXP (-ES/TS (I in*fcw*De/ *********************, 

USPAR (I) =u (NX, 0 ) ^ ' RS/ TSTAR*TS ( 1 , 1 ) /RSTAR/VEL 

CONTINUE 

RETURN 

END 
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NYP AR=NY + JT OP 

DO 120 1=1, 1 RIGHT 

begin CALCULATION 

DXl=DXP ( I ) , c * * * 

c*the ?Se 1 L owSo'S>s»m ? none, or state a» dxee. 

_ 1 °_^!iil^:I!----COMPUTE U-VELOCITY 

C "* DO 20 J=*l , NYPAR 

DYU= {DYP ( J+l) +DYP W \ \ ‘ 

DYL= (DYP ( J) +DYP ( J-l) ) / 2 • 

DY 2 =DYU+DYL rRVfl J) ) *DYL/DYU 

_DMOTV- ) -DTO/OTL ) /Dtt 

RHOVEE=DENS (1, j) * VPAR(1 ' Vr , XI 

ruoi=dens (l, J) *uparu, J) D 

RU02=RHOVEE*DYL/DYU/DY2 

rS°03=RHOVEE*DYU/DYL/DY2 

RU04=2.*GAV(1, J) /DY2/DYL 
RU05=2 . *GAV (1 , J) /DY2 /DYU 

RU06=DYL/DYU/DY2*DMUDY 

sr=s^-|f^™» 4 ;-»i:ssr BU01 

BB(J)= - RU02 03 ru04 -RU07 

CC ( J) ■“ _ T \ irR* f RINFNON— DENS (1 r J) ) _ 

2“l!ifl:-----SET BO^aS COITIONS FOR U-VELOCITY 

AA(0)=1- 

BB(0)=0. 

DD ( 0 ) =USPAR ( IBEG+I ) 

CC (NYPAR+1) =0 . 

CALL TDMA (0 , NYPAR+1 ) 

DO 30 J=0, NYPAR+1 

^——————————DETERMINE COEFFICIENTS FOR UPWIND SCHEME 

DO 3 9 U~ 1 f ^ \ then 

IF (VPAR(1, J+D -Vt.0.) THEN 

ALPH(J) =1- 

ELSE 

alph(J)=o. 

THEN 

BETA(J) =1 • 

ELSE 



142 


beta(j) = o. 

ENDIF 

39 CONTINUE 

~nn~In~T~7 COMPUTE T— 

DO 40 J=1 , NYPAR 

DYU=(Dyp (J+ i )+DYp(J))/2 
* DYP < J ) +DYP (J-l) ) /2 

DY2-DYU+DYL ’ 


_dmudy- ( + («*(!. J + 1, Mar ( i, j, , . DYW dyo 

RHOVEE-DEHS J)" 1 ’ ’ * D ™ /DrL > /°Y2 

SpSSOSSgSSr 

SawKs: &“? «& ' 2 • /di2 

RU06-DYL/DYU/DY2*DMubY 

RU07=DYU/DYL/DY2*DMUDY 

WFST— Q*DA*TPAR(1, J, 2) *TPARfl 7 •?» 

cc M {: +r „ 04 +ru05+ru » 6 /wMa.j.u 

40 DD(J)=RU0l*TPARn T -RU07 

C U,J,1)+WF ST*(3.-E/TPARn t n* 

AA (0 ) Zl SET B0UNDARY CONDITIONS FOR T--l~i 

BB(0)=0. 

DD (0) =TS (IBEG+1,1) 

CC (NYPAR+1) =0 . 

AA (NYPAR+1 ) =1 . 

DD (NYPAR+1 )=TPAR(1, NYPAR+1 l) 

CALL TDMA ( 0 , NYPAR+1 ) ’ 

DO 50 J=0, NYPAR+1 

50 TPAR?2 ( ^\'?V 0,) XX(J )=10- 

5 “ '?& ll <» W) . 1 - . a « J( . CT . ! . , 

_ DO 51 <J— 0 , NYPAR+1 

51 Sjyjl'aj ^/^Sr 0F STATE ^ DIFF - COEFF.'S*** 

EITOIf' =TPAR (2, J, 1) /TSTAR 

do _ 60~ j=o~nypar+i CUI<A ^ IE NEW VALUES of density 

c 60 DENS (2, J ) =TSTAR/TPAR(2, J, 1) 

~ T “ COMPUTE 

DO 70 J-l, NYPAR 

DYU— (DYP ( J+1 ) +DYP ( J) ) /2 

- MUDY ( + (GAT (l' J^GAT^ (1 ' J) ) *DYL/DYU 
RHOVEE=DENS d, jfivPARQ ’Z) 1 * * * DYU/DYL > /DY2 
RU01=DENS (1, J) *UPAR(1 wL 
RU22=RHOVEE*ALPH ( J) *2 ' /DY2 
RU23=RHOVEE* (BETA ( J) -Alph m i *, /1W „ 

RU2 4 =RHOVE E * (-BETA ( J) M2 /M2 * 7 


143 


RU 04 = 2 . *GAT ( 1 , J) / d * 2 /DYL 
RU 05 = 2 .*GAT( 1 , J) /DY 2 /DYU 

RU06=DYL/DYU/DY2*DMUDY 

rU 07 =DYU/DYL/DY 2 *DMUDY (2, J)*EXP (-E/TPAR (2, J, D ) 

S «“^™ 05 «U 06 -RU 07 -RU 0 8 

AA(J)-RU 01 + EU^ +RU 05 +RU 06 

g{?,: -RU 24 + RU 04 

_DD ,^-RU0i^ > Y . , . 

o****** ****** *******ASSUMED LEF LEO 

TE 1 =DENS ( 2 , 0 ) *VSPAR(IBEG + n*DYP(lW^ T (2#0|1) , 2 , TSTAR) ) 

_ ( GAMMAT(TPAR( 2 , 1 , 1 ) , 2 ,Ti>iAKj 

AA( 0 )= 1 .+TE 1 
BB( 0 )= 1 . 

DD ( 0 ) =TE 1 
CC (NYPAR+ 1 ) = 0 . 

DD (NYPAR+ 1 ) =TPAR ( 1 , NYPAR+ 1 , 2 ) 

CALL TDMA( 0 , NYPAR+ 1 ) 

DO 80 J ~ 0 f NYPAR+ 1 . / t\ pt 0 ) 

8 Q ?»R, 2 . 7 . 2 , -C^ S 7 ,. 0 ..XX^.O 1 ll 

DYL= (DYP ( J) +DYP (J-l ) ) / 2 • 

DY 2 =DYU+DYL CATfl J))*DYL/DYU 

_ DMUDI ' > -DW/D 7 I. I /D 72 

RHOVEE=DENS ( 1 , J) *vpar( 1 , J) 

RU 01 =DENS(l,J)*UPARa,J)/D X1 
PTl 22 =RHOVEE*ALPH ( J) *2 . /DYZ 9 

RU 23 =RHOVEEMBETA(J)-ALPH(Jn •/ 

RU 24 =RHOVEE* (-BETA(J) ) 2 ./DY^ 

RU 04 = 2 . *GAT ( 1 , J) /DY 2 /DYL 
RU 05 = 2 . *GAT ( 1 , J) /DY 2 /DYU 

RU06=DYL/DYU/DY2*DMUDY 

RU 07 =DYU/DYL/DY 2 *DMUDY * DENS ( 2 , J) *EXP (-E/TPAR ( 2 , J, 1 ) 

srsff i»SA!-“ 7 -r “ 8 

+ru 0 4 * RU ° 5 +RU 06 -RU 07 
CC(J)= -RU 24 +RU 04 

conditiors^or 

12 ■ /GWWAT ITPAR ' 2 ' 1 ' 1 ’ ' 3 ' TSTftR> 

AA( 0 )= 1.+ te 1 
BB ( 0 ) = 1 . 

DD ( 0 ) = 0 . 

CC (NYPAR+ 1 )= 0 . 

DD (NYPAR+ 1 > =TPAR ( 1 , NYPAR+ 1 , 3 ) 

CALL TDMA( 0 , NYPAR+ 1 ) 

S 

1 o 0 T? AR( 2 ,JjJ _ OOU UTE V-VELOCITY 

C VPAR ( 2 , 0 ) =VSPAR (IBEG+I) 
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DO 110 J=0, NYPAR 
DYU= (DYP (J+l) +DYP (j) ) /2 




IF 

IF 


(IBC.NE 


-SET B.C. 
1 ) GO TO 


S FOR 
114 


ELLIPTIC REGION- 


113 


DO 113 J=l, NY -° R ' { {IBEG - E Q-NX+1) .AND. (I.EQ.l) 

GVAVE= (GAV (1, J) +QAV (2 j) ) /•> 

. G ™= (GAT (1, J) +GAT (2 j )/2* 

Sc 3 REgfHi -£ AR (1 ' J) > /Dxi 

TBC ( J, 1) =GTAVE^TPAR (2 ???*.£'' J) > ^X1 
TBC (J, 2) =GTAVE* (TPAR (2 , j' ?l tdidm' 11 ' 1 * * /DX1 
TBC ( J, 3, =GTAVE* (TPAR a'. J.' 3) -TPM (J 5,' |J J 


) ) THEN 


ENDIF 


114 


TS ( IBEG+ 1,2) =TPAR^2 , If if 

KRITE 

WRITE (31, 900) 

WRITE (32, 900) 

WRITE (33, 900) 

WRITE (34, 900) 

WRITE (35, 900) 

IF (IDA.EQ.l) 


VALUES 


DATA 


112 


/™^ AR (2 , J) , J=0,NY+1+JTOP) 

TPAR (2 ' ? 1 ?"!' ^ +1+ JTOp ) 

TPAR 2 ' J'? ' ^bNY+I+jtop) 

THEN 2 ' J ' 3 ' J=0 ' NY+1+ JTOP) 

Dj E n2 A T-0 L w E r RSTAR/ASTAR 

U 112 J“ 0 , NY+I+Jrop 

WFPAR(J)=pref*den S(2 , J) *DENS(2,J)*TPAR(2 j ?1 
EXPf-E/TPAR/? t 1 \\ ' 


TPAR (2, J, 3) * 


CONTINUE 
WRITE (36, 900) 

ENDIF 

ENDIF 


(WFPAR ( J) , j=o, NY+1+JTOP) 


116 

120 

900 


no i,r 7 „ RESET VALUES 

DO 116 J=0, NY+1 + JTOP 

UPAR (1, J) =UPAR (2, J) 

VPAR ( 1 , J) =VPAR (2 , J) 

TPAR(1, J, 1)=TPAR( 2 , j i) 

TPAR (1, J, 2) =TPAR (2, J, 2) 

TPAR(1, J,3)=TPAR(2, j ,3) 

CONTINUE ' ' 1 

FORMAT (IP, 5 (IX, E12 5)) 

RETURN 

END 
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